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Abstract 



We derive the gravitational waveform and gravitational-wave energy flux 
generated by a binary star system of compact objects (neutron stars or 
black holes), accurate through second post-Newtonian order (0[(v/c) 4 ] ~ 
{S) ' 0[(Gm/rc 2 ) 2 ]) beyond the lowest-order quadrupole approximation. We cast 

the Einstein equations into the form of a flat-spacetime wave equation to- 
qq . gether with a harmonic gauge condition, and solve it formally as a retarded 

integral over the past null cone of the chosen field point. The part of this 
| integral that involves the matter sources and the near-zone gravitational field 

is evaluated in terms of multipole moments using standard techniques; the 
remainder of the retarded integral, extending over the radiation zone, is eval- 
*h . uated in a novel way. The result is a manifestly convergent and finite pro- 

cedure for calculating gravitational radiation to arbitrary orders in a post- 
Newtonian expansion. Through second post-Newtonian order, the radiation 
^ ' is also shown to propagate toward the observer along true null rays of the 

■ asymptotically Schwarzschild spacetime, despite having been derived using 

flat spacetime wave equations. The method cures defects that plagued previ- 
ous "brute-force" slow-motion approaches to the generation of gravitational 
radiation, and yields results that agree perfectly with those recently obtained 
by a mixed post-Minkowskian post-Newtonian method. We display explicit 
formulae for the gravitational waveform and the energy flux for two-body 
systems, both in arbitrary orbits and in circular orbits. In an appendix, we 
extend the formalism to bodies with finite spatial extent, and derive the spin 
corrections to the waveform and energy loss. 
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I. INTRODUCTION 



The generation of gravitational radiation is a long-standing problem that dates back to 
the first years following the publication of general relativity (GR). In 1916 Einstein calculated 
the gravitational radiation emitted by a laboratory-scale object using the linearized version 
of GR |l[ . Some of his assumptions were questionable and his answer for the energy flux was 
off by a factor of two (an error pointed out by Eddington Q]). There followed a lengthy debate 
about whether gravitational waves are real or an artifact of general coordinate invariance, the 
former interpretation being confirmed by the rigorous, coordinate free theorems of Bondi and 
his school [§],f||| and by the short-wave analysis of Isaacson |J . Shortly after the discovery of 
the binary pulsar PSR 1913+16 in 1974, questions were raised about the foundations of the 
"quadrupole formula" for gravitational radiation damping M (and in some quarters, even 
about its quantitative validity ||). These questions were answered in part by theoretical 
work designed to shore up the foundations of the quadrupole approximation P, [IT||n| , |T2| , |T3' 



and in part (perhaps mostly) by the agreement between the predictions of the quadrupole 
formula and the observed rate of damping of the pulsar's orbit fl^ . |l~5"|j . 

Because it is a slow-motion system (v/c ~ 10 -3 ), the binary pulsar is sensitive only to the 
lowest-order effects of gravitational radiation as predicted by the quadrupole formula. Nev- 
ertheless, the first correction terms of order v/c and {v/c) 2 to the quadrupole formula, were 



calculated as early as 1976 fi6|]17|| . These are now conventionally called "post-Newtonian" 
(PN) corrections, with each power of v/c corresponding to half a post-Newtonian order 
(1/2PN), in analogy with post-Newtonian corrections to the Newtonian equations of mo- 
tion In 1976, the post-Newtonian corrections were of purely academic, rather than 
observational interest. 

Recently, however, the issue of higher post-Newtonian corrections in the theory of gravi- 
tational waves has taken on some urgency. The reason is the construction of kilometer-scale, 
laser interferometric gravitational-wave observatories in the U.S. (LIGO project) and Eu- 
rope (VIRGO project), with gravitational- wave searches scheduled to commence around 
2000 (see |19| for a review). These broad-band antennae will have the capability of detect- 
ing and measuring the gravitational waveforms from astronomical sources in a frequency 
band between about 10 Hz (the seismic noise cutoff) and 500 Hz (the photon counting noise 
cutoff), with a maximum sensitivity to strain at around 100 Hz of Al/l ~ 10~ 22 (rms). 
The most promising source for detection and study of the gravitational-wave signal is the 
"inspiralling compact binary" - a binary system of neutron stars or black holes (or one of 
each) in the final minutes of a death dance leading to a violent merger. Such is the fate, for 
example of the Hulse- Taylor binary pulsar PSR 1913+16 in about 300 M years. Given the 
expected sensitivity of the "advanced LIGO" (around 2001), which could see such sources 
out to hundreds of megaparsecs, it has been estimated that from 3 to 100 annual inspiral 



events could be detectable f|l9lEU|,PT 



The urgency derives from the realization |22| that extremely accurate theoretical pre- 
dictions for the orbital evolution, and to a lesser extent, the gravitational waveform, will 
play a central role in the data analysis from these observatories. That data analysis is likely 
to involve some form of matched filtering of the noisy detector output against an ensemble 
of theoretical "template" waveforms which depend on the intrinsic parameters of the inspi- 
ralling binary, such as the component masses, spins, and so on, and on its inspiral evolution. 
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How accurate must a template be in order to "match" the waveform from a given source 
(where by a match we mean maximizing the signal-to-noise ratio)? In the total accumulated 
phase of the wave detected in the sensitive bandwidth, the template must match the signal 
to a fraction of a cycle. For two inspiralling neutron stars, around 16,000 cycles should be 
detected; this implies a phasing accuracy of 10~ 5 or better. Since v/c ~ 1/10 during the 
late inspiral, this means that correction terms in the phasing at the level of (v/c) 5 or higher 
are needed. More formal analyses confirm this intuition p3] , P^ , |2o] , [2B[ ] . 

The bottom line is that theorists have been challenged to derive the gravitational wave- 
form and the resulting radiation back-reaction on the orbit phasing at least to 2PN, or 
second post-Newtonian order, 0[(t>/c) 4 ], beyond the quadrupole approximation, and prob- 
ably to 3PN order. Furthermore, because of the extreme complexity of the calculations at 
such high PN order, independent calculations are called for, in order to inspire confidence 
in the final formulae. After all, the formulae will ultimately be compared against real data. 

This challenge was recently taken up by two teams of workers, one composed of Blanchet, 
Damour and Iyer (BDI), the other composed of the present authors. The goal was to derive 
the gravitational waveform and the energy flux for inspiralling compact binaries of arbitrary 
masses, through 2PN order. Each team adopted a different approach to the calculation, and 
worked in isolation from the other. Only at the end of the calculation were comparisons 
made for the key formulae for the waveform and the gravitational energy flux. The results 
agreed precisely |27| . 

The BDI approach was based on a mixed post-Newtonian and "post-Minkowskian" 
framework for solving Einstein's equations approximately, developed in a long series of pa- 



pers by Damour and colleagues p8| , p9| , |30| , |3T| , |32| , |33| . The idea is to solve the vacuum Einstein 
equations in the exterior of the material sources extending out to the radiation zone in an 
expansion ("post-Minkowskian") in "nonlinearity" (effectively an expansion in powers of 
Newton's constant G), and to express the asymptotic solutions in terms of a set of formal, 
time-dependent, symmetric and trace-free (STF) multipole moments [33]. Then, in a near 
zone within one characteristic wavelength of the radiation, the equations including the ma- 
terial source are solved in a slow-motion approximation (expansion in powers of 1/c) that 
yields a set of STF source multipole moments expressed as integrals over the "effective" 
source, including both matter and gravitational field contributions. The solutions involving 
the two sets of moments are then matched in an intermediate zone, resulting in a connection 
between the formal radiative moments and the source moments. The matching also provides 
a natural way, using analytic continuation, to regularize integrals involving the non-compact 
contributions of gravitational stress-energy, that might otherwise be divergent. 

The approach of this paper is based on a framework developed by Epstein and Wagoner 
(EW) [[RJ. Like the BDI approach, it involves rewriting the Einstein equations in their 
"relaxed" form, namely as an inhomogeneous, flat-spacetime wave equation for a field h a @, 
whose source consists of both the material stress-energy, and a "gravitational stress-energy" 
made up of all the terms non-linear in h a ^ . The wave equation is accompanied by a harmonic 
or deDonder gauge condition on h a/3 , which serves to specify a coordinate system, and also 
imposes equations of motion on the sources. Unlike the BDI approach, a single formal 
solution is written down, valid everywhere in spacetime. This formal solution, based on the 
flat-spacetime retarded Green function, is a retarded integral equation for h af3 , which is then 
iterated in a slow-motion (v/c < 1), weak-field (||^ a/3 || < 1 ) approximation, that is very 
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similar to the corresponding procedure in electromagnetism. However, because the integrand 
of this retarded integral is not compact by virture of the non-linear field contributions, 
the original EW formalism quickly runs up against integrals that are not well defined, or 
worse, are divergent. Although at the lowest quadrupole and first few PN orders, various 
arguments can be given to justify sweeping such problems under the rug JET), they are not 
very rigorous, and provide no guarantee that the divergences do not become insurmountable 
at higher orders. As a consequence, despite efforts to cure the problem, the EW formalism 
fell into some disfavor as a route to higher orders, although an extension to 3/2PN order 
was accomplished P5fl . 

One contribution of this paper is a resolution of this problem. The resolution involves 
taking literally the statement that the solution is a retarded integral, i.e. an integral over 
the entire past null cone of the field point. To be sure, that part of the integral that extends 
over the intersection between the past null cone and the material source and the near zone 
is still approximated as usual by a slow-motion expansion involving spatial integrals of 
moments of the source, including the non-compact gravitational contributions, just as in 
the BDI framework. But instead of cavalierly extending the spatial integrals to infinity as 
was implicit in the original EW framework, and risking undefined or divergent integrals, we 
terminate the integrals at the boundary of the near zone, chosen to be at a radius 1Z given 
roughly by one wavelength of the gravitational radiation. For the integral over the rest of 
the past null cone exterior to the near zone ("radiation zone"), we do not make a slow- 
motion expansion, instead we use a coordinate transformation to convert the integral into a 
convenient, easy-to-calculate form, that is manifestly convergent, subject only to reasonable 
assumptions about the past behavior of the source. This transformation was suggested by 
our earlier work on a non-linear gravitational-wave phenomenon called the Christodoulou 
memory |31S ]. Not only are all integrations now explicitly finite and convergent, we show 



explicitly that all contributions from the near-zone spatial integrals that grow with 1Z (and 
that would have diverged had we let TZ — > oo) are actually cancelledby corresponding terms 
from the radiation-zone integrals. Thus the procedure, as expected, has no dependence on 
the artificially chosen boundary radius TZ of the near-zone. In addition, the method can be 
carried to higher orders in a straightforward, albeit very tedious manner. The result is a 
manifestly finite, well-defined procedure for calculating gravitational radiation to high, and 
we suspect all, PN orders. 

The result of the calculation is an explicit formula for the gravitational waveform for a 
two-body system, the transverse-traceless (TT) part of the radiation- zone field, denoted h^, 
and representing the deviation of the metric from flat spacetime. In terms of an expansion 
beyond the quadrupole formula, it has the schematic form, 

h ij = {Q ij [l + O(e^) + 0(e) + 0(e^ 2 ) + 0(e 2 ) . . .j}^ , (1.1) 

where \x is the reduced mass, and represents two time derivatives of the mass quadrupole 
moment tensor (the series actually contains multipole orders beyond quadrupole). The TT 
projection operation is described below. The expansion parameter e is related to the orbital 
variables by e ~ Gm/rc 2 ~ (v/c) 2 , where r is the distance between the bodies, v is the 
relative velocity, and m = mi + wi2 is the total mass. The 1/2PN and 1PN terms were 
derived in ||17|1 , the 3/2PN terms in f35fl . The contribution of gravitational- wave "tails", 
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caused by backscatter of the outgoing radiation off the background spacetime curvature, at 
0(e 3 / 2 ), were derived and studied in ]32|j7] , 38 . 



This paper derives the 2PN terms including 2PN tail contributions; the results are in 

We also find that part of the tail terms at 3/2PN and 2PN 



complete agreement with BDI 39 



order serve to guarantee that the outgoing radiation propagates along true null directions 
of the asymptotic curved spacetime, despite the use of flat spacetime wave equations in the 
solution. The explicit formula for the general two-body waveform is given below in Eqs. 
( [STOP and (|6UTD . 

There are also contributions to the waveform due to intrinsic spin of the bodies, which 
occur at 0(e 3 / 2 ) (spin-orbit) and 0(e 2 ) (spin-spin); these have been calculated elsewhere 
|40|,|4l| , and are rederived in the EW framework in Appendix F. 

Equations of motion for the material sources must also be specified to 2PN order in order 
to have a consistent solution of Einstein's equations. These have the schematic form 



d 2 x/dt 2 = -(Gmx/r 3 )[l + 0(e) + 0(e 3/2 ) + 0(e 2 ) + 



[1.2) 



where x = x x — x 2 is the separation vector. The lowest-order contribution is obviously 
Newtonian. The next term 0(e) is the first post-Newtonian correction, which gives rise to 
such effects as the advance of the periastron. The term 0(e 3//2 ) comes solely from the spin- 
orbit interaction. The term of 0(e 2 ) is a second post-Newtonian correction to the equation 



of motion (and also contains spin-spin interactions). The terms in Eq. ( |1.2| ) are all non- 
dissipative, having nothing to do with gravitational radiation reaction. Through 2PN order, 
these equations are by now standard; see for example [|2|,[|3|,[|4[] and Eq. ( |6.5|) below. 



Given the gravitational waveform, we can compute the rate energy is carried off by the 
radiation (schematically J hhdQ, the gravitational analog of the Poynting flux). The result 
has the schematic form 



dE/dt = (dE/dt) Q [l + 0(e) + 0(e 3/2 ) + 0(e 2 ) + 



;i.s) 



Here (dE/dt)Q denotes the lowest-order quadrupole contribution, proportional to the square 
of three time derivatives of the trace-free mass quadrupole moment tensor of the source. The 
explicit formula for a general two-body system is given below in Eqs. ( |6.12|) and ( 6.13 ). For 
the special case of non-spinning bodies moving on quasi-circular orbits (i.e. circular apart 
from a slow inspiral), the energy flux has the form 



dE 
~dt 



2 /Gm\ 



32G 

'be 1 

/Gm\ 3/2 
+4W-rj + 



Gm Z2927 



re. 
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\ rc z J 



( 293383 380 
V 9072 + ~~9~ V 



:i-4) 



where i] = m 1 m 2 /m 2 . The first term is the quadrupole contribution, the second term is 
the 1PN contributon [|17]], the third term, with the coefficient 47r, is the "tail" contribution 
32l,|37lj38|,H5[l , and the fourth term is the 2PN contribution derived here. This new contri- 



bution was reported in [27], and was also derived using the BDI approach in |39]. For the 
contributions of spin-orbit and spin-spin coupling see PD| , |4T|J27| and Appendix F. 

Similar expressions can be derived for the loss of angular momentum and linear mo- 
mentum. These losses react back on the orbit to circularize it and cause it to inspiral. The 
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result is that the orbital phase (and consequently the gravitational-wave phase) evolves non- 
linearly with time. It is the sensitivity of the broad-band LIGO and VIRGO-type detectors 
to phase that makes the higher-order contributions to dE/dt so observationally relevant. 
For example, for an inspiral of two 1.4M Q neutron stars, the 2PN term in Eq. ( |1.4|) con- 
tributes about 9 of the 16,000 cycles observable in the bandwidth of the advanced LIGO. 
More detailed analyses of the effect of the 2PN terms on the matched filtering can be found 
in [25|,?fjf|7|. A ready-to- use set of formulae for the 2PN gravitational waveform template, 
including the non-linear evolution of the gravitational-wave frequency (not including spin 
effects) may be found in [Q. Spin corrections to the waveform templates may be found in 
Appendix F. 

An alternative approach to deriving gravitational waveforms and energy flux for inspi- 
ralling compact binaries, in the limit in which one mass is much smaller than the other, 
is that of black-hole perturbation theory. This method provides numerical results that are 
exact in v/c, as well as analytical results expressed as series in powers of v/c, both for non- 
rotating and for rotating black holes . For non-rotating holes, the analytical 



expansions have been carried to fourth PN order ||52|| . In all cases of overlap, the results 
agree precisely with our post-Newtonian results, in the limit rj —>■ 0. 

This paper is an attempt to present, in a relatively complete and self-contained form, 
the formalism and machinery of our "improved EW" approach to higher-order gravitational 
radiation from binary systems. Indeed, we begin with the raw Einstein equations, and end 
with a plot of the 2PN waveform. The goal is to provide sufficient detail to allow the reader, 
using this paper virtually alone, to verify any of the results reported here (we make no 
statement about the amount of work involved), and to carry the computations to higher PN 
orders. In Section II, we lay out the foundations of gravitational-wave generation, describing 
the relaxed Einstein equations, the matter sources and the near and radiation zones, and 
the formal retarded integral solution of the wave equation, including the new treatment of 
integration over the null-cone in the radiation zone. We turn in Section III to the weak-field, 
slow-motion approximation, and write down the matter and field variables to the accuracy 
needed to find the radiation to 2PN order. The part of the retarded integral for h af3 that 
extends over the near zone can be written in terms of a set of "Epstein- Wagoner" moments; 
these are evaluated explicitly in Section IV. In Section V, we evaluate the contributions 
to h a @ from the radiation-zone integrals, showing both the explicit cancellation of those 
terms in the EW moments that grow with TZ, and the generation of tail terms. Section 
VI specializes to two-body systems, and displays the full formulae for the gravitational 
waveform and energy loss. In Section VII, we further specialize to circular orbits. Section 
VIII makes concluding remarks. A number of technical details are relegated to Appendices. 

Our conventions and notation generally follow those of |53|,3~3}|. Henceforth we use units 
in which G — c — 1. Greek indices run over four spacetime values 0, 1, 2, 3, while Latin 
indices run over three spatial values 1, 2, 3; commas denote partial derivatives with respect to 
a chosen coordinate system, while semicolons denote covariant derivatives; repeated indices 
are summed over; rf v = rj^ = diag(— 1, 1, 1, 1); g = det(g^ u ); = (a*- 5 + a J *)/2; = 
(a u — a Jl )/2; e lJ is the totally antisymmetric Levi-Civita symbol (e 123 = +1). We use a multi- 
index notation for products of vector components: , with a capital letter 
superscript denoting a product of that dimensionality: x L = x n x n ...x 1 '; angular brackets 
around indices denote STF products (see Appendix A for definitions). Spatial indices are 
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freely raised and lowered with 5^ and 5ij. 



II. FOUNDATIONS OF GRAVITATIONAL- WAVE GENERATION 



A. The relaxed Einstein equations 

We begin our development of the gravitational- wave generation problem with the Einstein 
Equations 

Raf 3 _ _ g ^ R = 8vrT «/3 . (2.1) 

Here R a/3 is the Ricci curvature tensor, g al3 is the spacetime metric and T a/3 is the stress- 
energy tensor of the matter. Although Eq. ( |2.ip is a conceptually powerful statement, 
relating the curvature of spacetime on the left-hand side to the stress-energy of matter 
on the right-hand side, it is not a particularly useful form of the Einstein equations for 
practical calculations of gravitational- wave generation. For that purpose it is conventional 
first to define the potential 

h°0 = v a P - (-g) l/2 g aP , (2.2) 

(see e.g. |34j]) and to choose a particular coordinate system defined by the deDonder or 
harmonic gauge condition 

^%=0. (2.3) 

The spatial components of h al3 evaluated far from the source comprise the gravitational 
waveform and are directly related to the signal which a gravitational- wave detector measures. 
With these definitions the Einstein equations fl2.1|) can be recast in the following form 

nh af> = -167TT a/3 , (2.4) 

where □ = —d 2 /dt 2 + V 2 is the flat-spacetime wave operator. The source on the right-hand 
side is given by the "effective" stress-energy pseudotensor 

T a P = (-g)T al3 + (1671-)-^ , (2.5) 

where A alS is the non-linear "field" contribution given by 

A af3 = l^(-g)tf L + {h a ^ v bT) , (2.6) 

and t°^ L is the "Landau-Lifshitz" pseudotensor, given by 

16n(-g)tf L = {g x »g» p h a \ u h^, p + \g^h x \ p W\ v - 2 9lxu g x H^\ p h^ , x 

+ l -{2g aX g^ - g aP g^)(2g vp g aT - g pa g VT )h VT , x h^ J . (2.7) 
By virtue of the gauge condition (|2.3|) , this source term satisfies the conservation law 
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r a % = , (2.8) 

which is equivalent to the equation of motion of the matter T al3 .p = 0. 

We emphasize that Eq. ( |2.4j ) is not an approximate, or weak-field, form of the Einstein 
equations; it is exact, and relies only on the assumption that spacetime can be covered by 
harmonic coordinates. 

The form of Eq. (|2.4|) is suggestive of the wave equation for the vector potential in 
electromagnetism. This analogy with E&M is at once helpful and deceptive. It is helpful in 
that it suggests how to proceed to solve the equation, i.e. use a retarded Green function, 
and an expansion in terms of radiative multipole moments. It further illustrates that, just 
as the current density in E&M is the source for the vector potential, here the stress-energy 
of the matter is a source of the gravitational potential. 

However there are several important differences between Eq. ( |2.4| ) and its electromagnetic 
counterpart. First, the "source" in Eq. (|2.4|) also contains a gravitational part that depends 
explicitly on h a ^ , the very quantity for which we are trying to solve. Second, unlike the 
E&M case where the source (the currents) has finite spatial extent (compact support), we 
can expect r a/3 , which depends on the fields h a ^, to have infinite spatial extent. Indeed the 
very outgoing radiation that we hope to detect, will, at some level of approximation, serve 
as a contribution to the source, thus generating an additional component of the radiation. 
However, we have found that, for the physical situations of interest, this latter, highly 
nonlinear effect, often referred to as the Christodoulou memory, is very weak and can be 



adequately approximated by the methods of this paper |36| . 

Another complication in Eq. ( |2.4|) is that the second derivative term h al3 h^ u in the 
source really "belongs" on the left-hand side with the other second derivative terms in the 
wave operator. Such a term in a differential equation modifies the propagation characteristics 
of the field from the flat-spacetime characteristics represented by the d'Alembertian operator. 
Physically this is a manifestation of the fact that the radiation propagates along null cones 
of the curved spacetime around the source, which deviate from the flat null cones of the 
harmonic coordinates. Nevertheless, the techniques to be presented here do recover the 
leading manifestations of this effect, commonly known as "tails", including modification of 
the phasing of the solutions from their initial dependance on flat space retarded time to true 
retarded time of the asymptotic Schwarzschild spacetime of the source. 



B. Source, near-zone and radiation-zone 

We consider a material source consisting of a collection of fluid balls (stars) whose size 
is typically small compared to their separations. The material will be modeled as perfect 
fluid, having stress-energy tensor 

T af3 = {p + p)u a U P + pg a/3 , (2.9) 

where p and p are the locally measured energy density and pressure, respectively, and u a is 
the four-velocity of an element of fluid. We shall assume that the bodies are sufficiently com- 
pact that we can ignore all intrinsic multipole moments of the bodies at quadrupole order 
and beyond. That is, we treat only the bodies' monopole (mass) moments (in an Appendix 
we treat the bodies' dipole (spin) moments). For inspiralling binaries of compact objects, the 
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effects of rotationally induced and tidally induced quadrupole and higher moments on the 
orbital evolution or gravitational radiation have been shown, in the case of binary neutron 
stars, to be negligible until the final coalescence stage, where the post-Newtonian approxi- 
mation breaks down anyway ||54j| . For spinning black holes, the effects are small, but can be 



non- negligible for sufficiently large spin ||55| . In the long run, such finite-size effects should 
(and can) be incorporated into our formalism. 

To treat the monopole part of the bodies' mass distributions, we approximate the stress- 
energy tensor as a distributional tensor representing "point" masses, given by 

TSonopole = E^(-ff) 4/2 K"» 3 ( x " x aW) , (2-10) 

A 

where itia is the gravitational mass of the A-th body, and u\ is the four-velocity of its center 
of mass, x. A {t). Formally, such a distributional stress-energy tensor is not valid in general 
relativity. On the other hand, it has been shown in a variety of post-Newtonian contexts to 
give results that are equivalent to treating the bodies as almost spherical fluid balls, defining 
a suitable approximate center of mass, and carrying out explicit integrals over the interiors of 
the balls. The resulting self-field and internal energy effects result in a renormalization of the 
mass of each body from a "bare" mass J A pd 3 x to the gravitational mass tua- Furthermore, 
all effects of the internal structure of the bodies are "effaced", so that all aspects of the 
motion and gravitational radiation are characterized by a single mass itia for each body 
35| for demonstration of this effacement in the waveform at 3/2PN order). This is a 
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manifestation of the Strong Equivalence Principle, which is satisfied by general relativity. All 
these complications, then, can be embodied in the distributional stress-energy tensor of Eq. 
( p,10| ), with the caveat that all infinite self-field effects that might result from the use of the 
delta-function source are to be discarded (self-field effects having already been renormalized 
into tua)- An alternative viewpoint takes the gravitational field in a zone surrounding each 
body in a coordinate system that momentarily comoves with the body and notes that it can 
be characterized by multipole moments that can be identified with the body's asymptotially 
measured mass and (if desired) higher multipole moments. The fields surrounding each 
body are then matched to an appropriate interbody gravitational field, with the equations 
of motion providing consistency conditions for such matching. Apart from tidal effects, 
the results depend only on the effective masses of the bodies, and all self-field effects are 
automatically accounted for (see |56],[57] for example, for detailed implementations of this 



approach in various situations). 

The effects of spins can be added to the framework in a straightforward way; these are 
reviewed in Appendix F. 

We consider the bodies to comprise a bound system of characteristic size S = 
maxji^} tab, where tab = I x a — x B |, with a center of mass chosen to be at the origin 
of coordinates, X = 0. The source zone then consists of the world tube T = {x a \r < 
S, — oo < t < oo}. 

The bodies are assumed to move with characteristic velocities va < 1, and for much of 
their evolution with va <S 1. The characteristic reduced wavelength of gravitational radia- 
tion, X = \/2nr^S/v = 7l serves to define the boundary of the near zone, defined to be the 
world tube V = {x a \r < 1Z, — oo < t < oo }. Within the near zone, the gravitational fields 
can be treated as almost instantaneous functions of the source variables, i. e. retardation can 
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be ignored or treated as a small perturbation of instantaneous solutions. For most of the 
evolution, up to the point where the post-Newtonian approximation breaks down, TZ S> S. 

The region exterior to the near zone is the radiation zone, r > TZ. In this zone, we 
evaluate the fully retarded solutions of Eq. (|2.4j) , and focus on the parts that fall off as r _1 . 

The formal solution to Eq. Q2.4| ) can be written down in terms of the retarded, flat-space 
Green function: 



x) = 4 T ^^-^ix-xiy ^ (2 . n) 



r^(t / ,x')5(t'-t + 


x-x'|) 




x — x' 





This represents an integration of T a ^/|x — x'| over the past harmonic null cone C emanating 
from the field point (t,x) (see Fig. 1). This past null cone intersects the world tube T> 
enclosing the near zone at the three-dimensional hypersurface M . Thus the integral of Eq. 
( f2.11| ) consists of two pieces, an integration over the hypersurface A/", and an integration over 
the rest of the past null cone C — N '. Each of these integrations will be treated differently. 
We will also treat slightly differently the two cases in which (a) the field point is outside 
the near zone (Fig. 1), and (b) the field point is within the near zone (Fig. 2). The former 
case will be relevant for calculating the gravitational-wave signal, while the latter will be 
important for calculating field contributions to r a/3 that must be integrated over the near 
zone, as well as for calculating fields that enter the equations of motion. 

C. Radiation-zone field point, near-zone integration 

For a field point in the radiation zone, and integration over the near zone, we first carry 
out the t' integration in Eq. (|2.11|) , to obtain 







x — 


x' 


,x') 




X 


-x' 





feff(t,x)=4/ / V ; <fV • ( 2 - 12 ) 



Within the near zone, the spatial integration variable x' satisfies |x'| < TZ < r, where the 
distance to the field point r = |x|. We now expand the aZ-dependence in the integrand in 
powers of |x'|/r, using the fact that 



Ix-xT 



1 



E ^(-xV- lm (r q ),n.., m • (2-13) 



m=0 m[ 



We next expand T a/3 in a Taylor series about the retarded time u = t — r. The integration is 
now over the hypersurface A4, which is the intersection of the near- zone world-tube with the 
constant-time hypersurface t_M = u = t — r (see Fig. 3). Roughly speaking, each term in the 
Taylor series is smaller than its predecessor by a factor of order v < 1, thus for any hope of 
convergence of the series, one must restrict attention to slow-motion sources. We now have 
an infinite series in x' (expansion of |x — x'| _1 ) multiplying a double infinite series (expansion 
of |x — x'| inside the Taylor expansion). Grouping terms with the same powers of x' and 
carrying out the appropriate combinatorics (including use of "Faa di Bruno's formula" |f)8|j), 
it is straightforward to show that 

hff(t,x) =4E ^— p- (-M a/3k ^ k ") , (2.14) 



ql \r 
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where 

M a P kl - k «(u) = [ r a P(u, x')x ,fcl . . . x' kq d 3 x' . (2.15) 



M 

This general expansion, both in powers of r _1 and in retarded-time derivatives of 
M al3kl '" kq {u) will prove useful in later integrations of field quantities over the far zone. 

However, for gravitational-wave detectors, we need only to focus on the spatial compo- 
nents of h at3 , and on the leading component in 1/R, where R is the distance to the detector. 
Using the fact that = —N l , where N = X./R denotes the detector direction, we obtain 

4 00 1 f) m r 

m,*) = R^.O^ J M THuy)(K-xTd 3 x> + 0(R-Z) . (2.16) 
Because of the conservation law Eq. (|2.8| ), r lJ satisfies the identities 

r « = i(T 00 xV), 00 + 2(r'M),, - ~(r^V), w , (2.17a) 

T ij x k = I( 2t °(vV - r ofc xV'), + ^(2r l(i x j) x k - r*'xV"),, . (2.17b) 
Using these identities in Eq. (|2.16|) generates the multipole expansion 

o j2 00 

=R^T, ■ ■■N k jf^- km {u) , (2.18) 
where the "Epstein- Wagoner" (EW) moments are given by 

Ie\Y = J M T m X % X^d 3 X + IeW(sut{) , (2.19a) 

If w = J M (2r°^x k - r 0fe xV)d 3 z + J|^ (surf) , (2.19b) 
p ^... km = * / T ij x k, . ,. x ^ d z x ( m > 2 ) , (2.19c) 

TDj. (XL J J\/\ 

where integrating the spatial derivative terms in Eqs. ( |2.17|) by parts generates surface 
integrals at the two-dimensional coordinate sphere of radius 1Z bounding the hypersurface 
Ai, denoted dA4, resulting in surface contributions to the first two EW moments given by 

{d/dt) 2 4 w , {) = £ (4r'M - (r kl x^), k )n 2 h l d 2 n , (2.20a) 

JdM 

{d/dt)I% k ~ = £ (2r'Mx fe - ^x'x^tid 2 ^ , (2.20b) 

JdM 

where h l denotes an outward radial unit vector, and d 2 Q denotes solid angle. 

One advantage of this multipole expansion is that the field and source variables appearing 
in the integrand r a/3 are evaluated at the single retarded time u\ a disadvantage is that 
because the field contributions to r a/3 fall off as some power of r, one can expect to encounter 
integrals that depend on positive powers of the radius 1Z of the boundary of integration, 
especially in some of the higher-order moments. If this boundary is formally taken to oo 
(as was previously done), these integrals would diverge. However, as we shall see, such 1Z- 
dependent effects are precisely cancelled by contributions from the integral over the rest of 
the past null cone, to which we now turn. 
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D. Radiation-zone field point, radiation-zone integration 



The integral over the rest of the past null cone C — Af can be written in the form 



hfjt, x) = 4 / du' 1 ' 1 > i ' — S(u' - 1' + r ') d ^' , (2.21) 

J-oo JC-N |X — X'| 

where we have simply inserted 1 = / du'8{u' — t' + r'). We now integrate over t' and r', and 
note that 

I" dt' r dr'5{u' -t' + r')5(t' - t + |x - x'l) = / t-u'-J- x u> < u and f > 11 ( 2 .22) 
J-oo Jn [0 u > u or r < 7c . 

The result is 

hgjr(t,x) = 4 f d«' / ^^±^0 [^,nO]W . (2.23) 
J-oo Jc-N t — u — n • x 

Note that r' is a function of v! and f2' via the condition [from the two delta-functions in Eq. 
fl2.22 )]: t — u' — r' + |x — x'|, which gives 

r\u', n') = [(t - m') 2 - r 2 ]/[2(t - u' - n' ■ x)] . (2.24) 

The integration over solid angle d 2 Q' for a given value of u', together with the v! + r' 
"time" dependence of r a/3 , can be seen to represent an integration over the two-dimensional 
intersection of the past null cone C with the future null cone t' = u' + r' emanating from 
the center of mass of the system at tcM — u' (Fig. 4). The integration over u' then includes 
all such future-directed cones, starting from the infinite past, and terminating in the one 
emanating from the center of mass at time u, which is tangent to the past null cone of the 
observation point. 

However, for u > u' > u — 2TZ, the two-dimensional intersections meet the boundary of 
the near zone, and so the angular integration is not complete. If we choose the field point 
x to be in the z-direction, so that n' ■ x = r cos 8', then the condition r' > TZ, together with 
Eq. (ggg ) imply that < <p' < 2vr, 1 - a < cos#' < 1, where 



a. 



(u-u')(2r-2K + u-u , )/2rK . (2.25) 



Note that a ranges from {v! — u) to 2 (v! = u — 27c). For u' < u — 27c, the angular 
integration covers the full 4tc. Thus we write the radiation-zone integral in the form 



i-u f'2m r L t" p In '4- r r, 

hf_ N (t,x) = A du' d(f>' { [' J [r'(u',n')] 2 dco86' 



u-2Tl Jo Jl-a t — U 1 — fl' • X 

+4 / du' I { , !; } [r'{u', Q')] 2 d 2 n' . (2.26) 



DC 



t — u' — n' • x 



Note that r a/3 contains only field contributions evaluated in the radiation zone; in determin- 
ing these we will make use of the general expansion ( |2.14j ). 

To obtain the contribution to the gravitational waveform, we evaluate the spatial com- 
ponents of Eq. Q2.26D at distance R and direction N and keep the leading 1/R part. 



12 



E. Near-zone field point, near-zone integration 



In this case, in Eq. ( |2.12| ), both x and x' are within the near zone, hence |x — x'| < 21Z. 
Consequently, the variation in retarded time can be treated as a small perturbation, since 
T af3 varies on a time scale ~ 7Z. We therefore expand the retardation in powers of |x — x'|, 
to obtain 

00 1 r) m r 

4 / r^(t,x')|x-xr-W, (2.27) 

where M. here denotes the intersection of the hypersurface t = constant with the near-zone 
world-tube. 



F. Near-zone field point, radiation-zone integration 



The formulae from Section IIP , such as ( (2.24|) and (p.25|) , carry over to this case with 
only one modification. The final future null cone that appears in the integration is the 
one that intersects the boundary of the near-zone and the past null cone of the field point 
simultaneously at u' = u — 27Z + 2r, rather than u' = u (Fig. 5) (recall that here, r < 71). 
The result is, for a near-zone field point, 



C-M 



(t,x) 



-21l+2r 



dv! I dip 



* r"V + r',x') r Jt . J o/M2. 



u-27L JO Jl-a t — U' — II' ■ X 

r u-2K f T <x$(o,l _|_ r * 



[r'{u',n')] 2 dcQs6' 



[2.28) 



G. Gravitational waveform and energy flux 

To obtain the gravitational waveform, we combine the two contributions to h % \ Eqs. 
( j2. 181 ) and the leading 1/R part of the spatial components of (|2.26| ), and evaluate the 
transverse-traceless (TT) part, given by 

Kt = h k \nP? - \p^P kl ) , (2.29) 

where P l k = 5 l k - N k N\ 

Note that the two expressions that contribute to h kl in Eq. ( [2.29P each depend on the 



radius 1Z of the near zone. Since 1Z was an arbitrarily chosen radius, the final physical answer 
should not depend on it. However, to check that all terms involving 1Z cancel in the end 
would be a formidable task. Instead we adopt the following non-rigorous, but reasonable 
strategy. All terms in the near- zone EW moments and in the radiation-zone integrals that are 
independent of 7Z are kept. All terms that fall off with 7Z will be dropped. Close examination 
shows that, despite our formal choice 7Z ~ A, nothing in our calculations actually constrains 
the value of 7Z, apart from the inequality 7Z < R. Thus we are free to make 7Z sufficiently 
large, but still less than R, so as to make such terms as small as we wish, whether or not 
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they ultimately cancel. In this regard, it is useful to note that, for a LIGO/VIRGO detector 
10 Mpc from a source emitting gravitational waves at / = 100 Hz, fR = R/X ~ 10 16 , and 
thus many orders of magnitude of 1Z are available to achieve this suppression. Nevertheless, 
we believe that all such terms actually cancel. Finally, all terms that grow with powers of 
TZ are also kept. In this case we will show explicitly that all terms that vary as positive 
powers of TZ cancel between the near-zone and radiation-zone integrals. This procedure thus 
isolates the finite terms that arise from convergent integrals, while simultaneously verifying 
that no truly divergent integrals arise. The result is a well-defined, explicitly finite, method 
for calculating the gravitational waveform. It is the explicit inclusion of the radiation-zone 
integral in the formulation of Eq. ( |2.26| ) that cures the apparent divergences that plagued 
the original EW framework. 
The energy flux is given by 

E = (R 2 /32tt) £ h^ T h% T d 2 Vt . (2.30) 



III. WEAK FIELD, SLOW-MOTION APPROXIMATION 

A. Iteration of relaxed Einstein equations 

We make the standard assumption that, with respect to the orbital motion and mutual 
gravitational interactions, 

v\ ~ ttia/S ~ e < 1 . (3.1) 

where e will be used as an expansion parameter. 

Now, because the field h alS appears in the source of the field equation, the usual method of 
solution is to iterate: substitute h af3 = in the right-hand side of Eq. ( 2.11|) and solve for the 



first-iterated hi ; substitute that into Eq. (|2.11|) and solve for the second-iterated h% , and 



so on (imposing the gauge condition Eq. ( |2.3| ) consistently at each order). The first iterated 
hi 13 is 0(e), and each subsequent iteration improves its accuracy by one order in e. Thus, 
for example, to obtain a result for the waveform accurate to the order of the quadrupole 
formula, h ~ (m/r)I^ ~ (m/r)(v 2 + m/ S) ~ e 2 , two iterations of Eq. (|2.11|) are needed. To 



obtain the first post-Newtonian corrections to the quadrupole approximation, i.e. h to order 
e 3 , h^f , or three iterations, are needed, while to obtain the 2PN contributions (the goal of 
this paper), the fourth-iterated field is needed. This would be a daunting task, if it weren't 
for the use of the identities, Eqs. ( |2.17| ). Consider for example, the quadrupole formula. The 



source of the second-iterated field h % l contains pv % v J as well as terms of the form (V/^ )" 
both of which are O(pxe). (Note that (Vh) 2 ~ KW 2 h ~ pe). However, the use of the identity 



Eq. ( |2.17a ) in the near- zone integration converts r*- 7 into two time derivatives of r 



00 fj^i % 



(modulo total divergences); because of the slow- motion approximation, two time derivatives 
increase the order by e, and thus, to sufficient accuracy, only the dominant contribution to 
r 00 , namely p, is needed, without explicit recourse to the first-iterated h^ . Instead, h^f is 
buried implicitly in the equation of motion ( [2.8|) that leads to the identity ( |2.17a|) . This 



circumstance is responsible for the prevalent, but erroneous view that linearized gravity (one 
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iteration) suffices to derive the quadrupole formula. The formula so derived turns out to be 
"correct", but its foundation is not (see || for discussion). 

Thus, in practice, in order to evaluate EW moments required for the N-th iterated field, 
we will only need the (N-2)-iterated field contributions to the sources. This is not precisely 
true for the two EW surface integrals, and formally the full (N-l)-iterated field must be used 
in t u there, but with sufficient care, it can be shown without detailed, explicit calculations 
that the contributions of the (N-l)-iterated fields all fall off sufficiently rapidly with 71 to 
have no effect on these surface integrals. Similarly, for the radiation-zone integration, the full 
(N-l)-iterated field must be used in r y , Eq. ( f2.2tj| ). However, it will also be possible to show 
that the contributions of the these fields fall off with 7Z. To obtain the finite contributions 
and the contributions needed to cancel any divergent terms from the EW moments, only 
the N-2 iterated fields will be needed in practice. Thus to 2PN order (fourth iteration), only 
second-iterated fields will be needed explicitly in the source terms. 



B. Second-iterated fields in source terms 

Because the source contributions are integrated over all space, we must evaluate the 
second-iterated fields h% in a form that is valid everywhere (this and the following section 
follow the approach and notation of BDI; see P3"|, for example). The first iteration of 



ah 00 = 


— 167T 


ah 0i = 


— 167T 


oh ij = 


— 167T 



the field equations (|2.4| ) gives the linearized equations, Oh" = —16iiT al3 . Since T a P has 
compact support, the solutions are standard Lienard-Wiechert-type retarded functions. The 
solutions have the leading order behavior h 00 ~ e, h° l ~ e 3//2 , h 1 ^ ~ e 2 . Taking these 
orders into account, together with the fact that, because of the slow-motion assumption, 
d/dt ~ e l l 2 d/dx % , we can write the second-iterated field equations in the form (we drop the 
subscripts) 

-g)T™+ 7 -h™h™ + 0( P e>), (3.2a) 
-g)T 0i + 0(pe 3 / 2 ) , (3.2b) 
-g)T« - \{h™h™ - U l3 h™h™) + 0( P e 2 ) , (3.2c) 

where we have kept only contributions required to determine h 00 , h 0t , and to the ac- 
curacies e 2 , e 3//2 , and e 2 , respectively (note that, in identifying orders of source terms with 
dimension (length)" 2 , we can use 0~ 1 p ~ e). By defining the densities 

o = T 00 + T u , (3.3a) 

^ = T 0i , (3.3b) 

a l3 = T ij , (3.3c) 

and the retarded potentials 

V(t,x)=[ , dX , cr(t — Ix-x^x 7 ) , (3.4a) 
Jc |x — x'| 

Vi(t,x)= I -^-^(t-lx-xVO, (3.4b) 
Jc x — x' 
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Wy(t,x)= / 



d 3 x' 



c x 



(t 



x — X 



it is straightforward to solve Eqs. ( |3.2|) to the needed order, with the result 



h 00 = AV - 4(W - 2V 2 
h 0i = AVi + 0(e 5/2 ) , 
h ij = 4W ij + 0(e 3 ) , 



+ 0(e 3 



(3.4c) 



(3.5a) 
(3.5b) 
(3.5c) 



where W = Wu. It is useful to note that, although these forms of h a ^ are of sufficient 
accuracy in practice to be used in the effective sources for evaluating the waveform to 2PN 
order, they are not sufficiently accurate for use in the equations of motion that must also 
be specified consistently to 2PN order. The 2PN equations of motion require h 00 to 0(e 3 ) 
and h° l to 0(e 5 / 2 ) (h^ is sufficiently accurate as it stands). However, as the 2PN equations 
of motion are well known, we shall not undertake their derivation here, and will simply use 



the published equations [@,^] when they are needed. 

Because the source of V and V{ has compact support, the integrals ( |3.4a| ) and ( |3.4b| ) can 
be evaluated simply for field points within either the near zone or the radiation zone. But 
because the source of Wij contains both compact and non-compact support pieces, it must 
be evaluated carefully, with proper attention paid to contributions from the integration over 
the radiation-zone part of the null cone. The details will depend on the use to which Wij is 
being put. Evaluation of is discussed in Appendix C. 

When we calculate the EW moments, we shall need the field contributions to r a/3 eval- 
uated at fixed retarded time u (on the hypersurface A4), and for field points with r < 1Z. 
We therefore expand the retardation t — |x — x'| as a perturbation of the potentials V, Vi 
and Wij about t = u, with |x — x'| acting as the expansion parameter [see Eq. (|2 .27|) ] . The 
results are 



V 



U + -d 2 X + 0(e 5 / 2 ) 



W l3 



Vi = Ui + 0(e 5 / 2 ) , 



where the "instantaneous" potentials are given by 



d 3 



-<7[U. X 



J M |X — X'| 

X(u,x)= / d 3 x'\x. — x.'\a(u, x') , 

J M 



d 3 



X' 



.M X — X' 



d 3 



X' 



M X — X' 



CTi[U, X 



-5ijU tk U tk ) 



{ u, X 



(3.6a) 

(3.6b) 
(3.6c) 

(3.7a) 
(3.7b) 
(3.7c) 
(3.7d) 



We have used the fact that, by virtue of the conservation of mass and momentum at lowest 
order, d t J ad 3 x ~ e 5//2 and d t J Oid 3 x ~ e 3 . We will drop the contribution from the radiation- 
zone integral (Wij)c-Ni which falls off at least as fast as 1Z~ 2 (see Appendix C). Note that 
these potentials satisfy U iti — —If, V 2 A = 2U, Pijj — —Ui. 
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C. Near zone metric, matter stress-energy, and effective gravitational source 



In order to evaluate the components of the stress-energy tensor T af3 to the necessary 
order, we need the components of the near-zone metric to post-Newtonian order. These are 
given from Eqs (|2.2|) and (|3.5|) by 



no 



'1 + 2V + 2V 2 ) + 0(e 3 



,0i 



-4^ + 0(e 5 / 2 ), 
1 - 2V)8 ij + 0(e 2 ) 



(3.8a) 
(3.8b) 
(3.8c) 

(- g ) = \ + 4V -8(W -V 2 ) + 0(e 3 ) . (3.8d) 

These equations, together with the distributional definition ( |2.10|) of the stress-energy tensor 
yield, to the requisite order, 

3 



a 



A 



+ l -V 2 + l -Vv\ + AW+ ? -v 4 A - AV t v\ + 0(e 3 ) 



5 3 (x-x A ) 



a 



A 



l-V + -v 2 A + 0(e 2 ) 



5 3 (x-x A ) , 



l-V+ l -v 2 A + 0{e 2 ) 



(3.9a) 
(3.9b) 

(3.9c) 



where the potentials V, Vi and W are assumed to be evaluated at x A , excluding contributions 
of the A-th body itself (to avoid infinite self-field terms). The components of T al3 can be 
easily constructed from these expressions. 
To the needed order, h. a " has the form 



A 00 = -UV k V k 



16 



-vv + v k v k - 2v k v k + |y 



1 7 

+-Kn,fe(Kn,fc + 3\4, m ) + 2W, k V k - W k iV k i - -vv k v k 



+ 0(pe 3 



A 



0/ 



16 



V k (V k . 



v Lk ) + Ivv. 



4 



+ 0(pe 5 / 2 ) 



A y = 4 ( ViVj 



SijV k V k ) + 16 



2V {i V j) -V k M 



+2T4,(iV^) )& — 8ij(—V + V k V k — V mik V[ mik ] 



k,i v k,j — Vi,kVj,k 

+ 0(pe 3 ), 



(3.10a) 
(3.10b) 

(3.10c) 



where overdot denotes d/dt. Notice the presence of the cubically nonlinear terms in A 00 , 
involving either V x W or V 3 . 



IV. EVALUATION OF EPSTEIN- WAGONER MOMENTS 

A. Basic strategy 

The EW moments are integrals over a sphere of harmonic coordinate radius TZ about 
the center of mass of the system, with all variables entering the integrands to be evaluated 
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at retarded time u = t — r. We substitute the matter stress-energy tensor T Q/3 , and the 
second-iterated fields evaluated in the near-zone into Eqs. ( 2.19Q . We expand all quantities 



to the PN order needed to achieve a 2PN-accurate waveform. Each volume integral will be 
split into a "compact" (C) piece involving integration of the compact-support matter source, 
and a "field" (F) piece, involving integration of the non-linear field contributions. In I^ w 
and the two surface integrations at the boundary radius 1Z will involve only the field 
contributions, and will require somewhat special treatment. 

In integrating the field terms, we will frequently integrate by parts, but will carefully 
evaluate and save the surface terms, using the identity 



M 



dkF v-m d s x = a F l] - m \ n n k Tl 2 d z n . (4.1) 

M JdM 

In order to simplify some of the integrations, we will frequently make a change of variables 
within integrals, in order to place one of the bodies at the origin of the new variables, for 
example y = x-x^. Even though d 3 y = d 3 x, this shift has the consequence that the region 
of integration M. x = {x l \\x.\ < 71} will now appear in the new coordinates to be a region 
bounded by |y| = \7Zh — x A |, i.e. not centered at y = 0. It is much easier in practice to 
integrate in ^/-coordinates over a region M. y = {y*||y| < 7Z}, which is shifted by x^ relative 
to the true region of integration. The two integrations can be related by taking into account 
the appropriate surface integrals, using the identity 

/(x)d 3 x = / g(y)d 3 y - <f g(y)y ■ x A K 2 d 2 Q v 

J My J8My 

+-l x A ■ V<?(y)y ■ ^ A 7Z 2 d 2 Q y + ... , (4.2) 

£ JdMy 

where g(y) = f(y + x^) and y = y/y. Again, we evaluate and save the surface terms. 

In the end, we will only be interested in the physically measurable, transverse-traceless 
(TT) components of the radiation- zone field h %3 . We will therefore make frequent use of the 
identities, which follow from the definition (|2.29| ): 

(S ij )tt = , (A^')tt = , (4.3) 

where B is arbitrary. These identities apply only to indices "i" and "j" appearing in the 
components of the final waveform; we do not apply them to fields which ultimately make 
up source terms. 

In the field integrals, we will need explicit forms for the instantaneous potentials ( p.7|) 
evaluated inside the near zone. To the needed order, they are given by 

(4.4a) 
(4.4b) 

(4.4c) 

i ^ r + ^ 3 ), (4-4d) 

x - xa\\x a - x B 



U(u, 


*) = 5 




- + 0(e 3 ), 




X l x - X A 


X(u, 


*) = 5 


^m A |x - 

4 


x A |(l + 0(e)), 


Ui(u, 




m A v l A 

X l x - X A 


- + 0(e 5 / 2 ), 


P(u, 


*) = 5 


- m A v 2 A 


4 2 ^ 


X l X - *A 
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where P 



Pa, and where 



m A = m A 




^m B /|x A - x B | + 0( 



(4.5) 



B 



Equation ( 4.4d|) can be easily obtained from Eq. ( |3.7d| ) (after contraction on ij) by inte- 
grating by parts, carefully checking the vanishing of all surface terms. Although the full 
potential Pjj appears (via Wij) in A 00 , we will not need its explicit form, as the integra- 
tion of that particular term will be handled by a "trick" (see Appendix D). Note that the 
so-called "superpotential" X(u, x) is needed only to lowest order because it always appears 
twice time-differentiated, e.g. in Eq. ( |3.6a| ), and so its contribution is already 0(e) relative 
to that of U. 



B. The two-index moment P^W 

We write Eq. ( |2.19aj) in the form 

1^ = 1^ + 1^ + 1^ (4.6) 

where the three terms represent the compact (C), field (F) and surface (S) contributions. 
Substituting Eqs. @, (TO , ( ggp , flSTSD , and flOj) into (-g)T 00 and expanding through 
0(pe 2 ), we obtain 

r£ = Y. m Axl ( 1 + \v\ + 3^2— ) + lj2 m A^A v4 



2 b r AB 



A 



x ij ( 7 3 

+ ^m A m B — 2v 2 B + -v\ - 4v A • v B - -(v B ■ rUs) 2 - 
AB r AB \ 2 2 c 



r BC 



+lY, — -l a B-^AB)+0(e 3 )xmx 2 A , (4.7) 
2 c r AC 2 J 

where n. AB = x^ — x B , r AB = \x.ab\, ™-ab = * AB /r AB , and = d 2 ^ A /dt 2 . All sums are 
assumed to exclude cases where a denominator (e.g. r B c) might vanish. 

To the required order for calculating Pp , A 00 can be written in terms of the instantaneous 
potentials, 



A"" I [U k U k + 16 {~U k X k -UU + U k U k - 2U k U k + ^U 2 

+\u m>k (U m , k + 3U k , m ) + 2P, k U k - P km U km - ~UU k U k } + 0(pe 3 ) . (4.8) 

For the first term, the integral — (14/1671") f M U jk U jk x t: >d 3 x is straighforward: integrating 
twice by parts and showing that the surface terms are proportional to 715^ , which has no TT 
part, we are left with the integral (14/167r) J m UV 2 Ux %:) d 3 x = — (7/2)i^ AB m A m B x l A /r AB . 
This term is of 1PN and 2PN order via the PN contributions to m*. The next term, 
— (14/1671") J M U jk X tk x^d^x is already of 2PN order. We integrate once by parts to re- 
move the derivative from U. Using the fact that V 2 X = 2U, we find a surface in- 
tegral -(14/1671-) § dM UX^ k x ij n k n 2 d 2 n, and the new integrals (28/16tt) J m UUx ij d 3 x + 
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(28/167r) J M U X ^8 k ^ % x^ d 3 x . The first of these volume integrals can be combined with that 
arising from the third term in Eq. ( |4.8|) . We substitute Eqs. ( |4.4a ) and ( |4.4b|) , including 
a delta-function term that arises in U (see Appendix B). In the surface term, we expand 
the integrand in powers of r _1 , and obtain —(7/15) Y,AB m A lr nB / TZWi + x^a*}) + 0(1Z~ 1 ). 
We drop all terms that fall off with increasing 1Z. In the volume integrals, for each term 
in the sums U = J2aUa and X^ = J^A^A,k, we change integration variables from x to 
y = x — x^ so that, for a given A, the potentials U a and X A ,k are centered at the origin 
of the new y coordinate, while U now takes the form J2B m B/\y + %ab\- We calculate the 
surface contributions that result from this change of variables using Eq. ( f4.2|) . For example, 
the first integral then becomes 



/ UUx lj d 3 x = ^2m A m B 



I My ly + x^ 

x(y 2 y^ + 2yy^2 + xfjy- 3 d 3 y. 

We use the spherical harmonic expansion 

1 



ya A ■ y - v\ + 3(v A ■ y) 2 - y^y 3 5 3 (y) 



l, 711 



(4.9) 



(4.10) 



where r<(>) denotes the lesser (greater) of tab and y, express all products of unit vectors 
y k in terms of symmetric, trace-free (STF) products using Eqs. (|A2|) , and integrate over 
directions y, using the identity 



E / Y; m {h)Y lm {y)y <v> d 2 n y 



n 



<L> , 



(4.11) 



(see Appendix A) where the superscript < L > over a unit vector denotes an /-dimensional 
STF product. We then integrate over y, using the formula 



o r 



i+i 



y q dy 



{2i + 1; 



AB 



(l + q + l)(l-q) 



1 - 



l+q+1 ( n 



21 + 1 \r A B 



q -i 



(4.12) 



The result is a series of terms of three types: those with non-vanishing TT part that are 
independent of 1Z and linear in 1Z, which we keep; terms with vanishing TT part which we 
discard (regardless of their dependence on 1Z); and terms that fall off with increasing 7Z, 
which we also discard. An example of the second type of term would be a contribution to P^ w 
proportional to 7Z5 lJ . The contribution of such a term to h l i has no TT part; equivalently, 
it can be eliminated to the necessary order by a finite gauge or coordinate transformation. 

Many of the field integrations that we encounter in evaluating the EW moments are 
amenable to this general method: (i) integrate by parts to leave one potential undifferen- 
tiated, (ii) change variables to put the center of the differentiated potentials at the origin, 
(iii) expand the undifferentiated potential in spherical harmonics, (iv) express all unit vec- 
tor products in STF terms, (v) integrate over d 3 y using the identites ( |4.11|) and (|4.12|) , (vi) 
retain all relevant contributions from surface integrals that arise in steps (i) and (ii). 

Terms 2 through 8 contributed by A 00 [Eq (|4.8p ] can be handled using this method, as 
can the compact contributions to Pjj and P (proportional to velocities) in terms 9 and 10. 
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However the non-linear field contributions to Pjj and P lead to additional complications, 
although the basic method still applies. These terms are discussed in Appendix D. Finally, 
the cubically non-linear term 11 can be calculated easily by integrating by parts. Compu- 
tation of these terms is straightforward but tedious. In evaluating 2PN terms, we make 
repeated use of the fact, valid to Newtonian order, that Y^A m AX-A = 0. 



We now turn to the surface term To, given by Eq. (|2.20a|). Because the surface lies 



L s > 

outside the matter source, only the field contribution, A*- 7 ' is needed. The term can be 
rewritten in the form 

(d/dt) 2 ii j = (i/i6tt) / hA^h^n 3 - A kl ^n 4 ) d 2 n . u.u) 

JdM V ' 

However, because Ig is essentially two anti-time-derivatives of the surface integral, re- 
ducing its order by e, we need to know A*- 7 to 0(pe 3 ), i.e. to 0(e 2 ) beyond its leading order 
terms, at least in principle. This is in contrast to having to know A 00 in the spatial integral 
Ip only to 0(e) beyond its leading order. This would present considerable complications, 
except for the fact that we only need to calculate a surface integral, and retain terms that 
are either independent of or grow with 1Z. Consequently we only need to retain contribu- 
tions to A'- 7 that vary as 1Z~ 2 or 1Z~ 3 . To see what terms must be retained, we return to 
the definition of A*- 7 , Eq. (|2.6|) . Far from the source, the fields h af3 have the leading e and 
r dependences h 00 ~ e/r, h° l ~ e 3 / 2 /r 2 (r~ 2 here because the net momentum of the system 
vanishes), and h l i ~ e 2 /r; A*- 7 has the schematic form {h\) 2 + h(h \) 2 + h 2 (h^\) 2 + . . .. By 
combining the leading forms of h a/3 with the knowledge that time-derivatives increase the 
order by e 1 / 2 , while spatial derivatives either increase the rate of fall-off by one power of r -1 
or increase the order by e 1 / 2 via the retarded time dependence, it can be shown by inspection 
that terms of order h(h^\) 2 and higher are either of higher than 2PN order, or fall off faster 
than TZ~ 3 , or generate angular dependence that leads to no TT parts. However, the purely 
quadratic terms proportional to (h^) 2 do contribute; their explicit contribution is given by 
the non-linear terms of Eq. (|2.6|) with g^ u replaced by 77^. Again, inspection shows that, to 
the required order, we can write 

A ii = _ h nftj + Ifcoo^ + 2fc 00 '(W> - ^(h 00 , k h 00 ik + h 00 > k h k0 ) . (4.14) 

Further inspection shows that knowing h a/3 to the accuracy shown in Eq. ( |3.5| ) suffices; 
the higher-order terms not explicitly shown in those expressions contribute terms either at 
higher-than-2PN order, or at faster-than-7?. -3 fall-off. We do need to evaluate V, Vi and 
Wij carefully, however. Expanding these functions in powers of |x — x'| about t = u, but to 
higher orders than that shown in Eq. ( j3.6|) , using Eqs. ( |3.9|) for er, <jj and a^, and displaying 
only terms that lead to the appropriate contributions in A u , we find in the vicinity of r = 1Z 

V = ^ + ±Q kl (36 kl - n k n l ) + 0(e 3 /r) - ~ Q +0(e 3 r°) 
r 4r 3 



, (4) 

+ — Q kl {S kl + n k n l ) + 0(e 4 r) + 0(e 2 /r 2 ) + 0(e 2 /r 3 ) , (4.15a) 

(3) 



V t = — -U ljk J k - Q lj W + 0(e 5/2 )/r 2 - 4 Q ij ri + 0(e 3 r°) + 0(e 3/2 /r 3 ) , (4.15b) 
2r 2 4 
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X 


= to 


1 ^ m^XA 1 + 

A V 


Q 11 


= E 

A 




Q ljk 


= E 


m A xi , 


J' 


= E 


m A t ^a v A i 




= E 


m A e dm x l A v A x j A , 



W ij = }_Qij + {e 2 /r 2 ) + 0(e 3 /r) , (4.15c) 

where we define here and for future use 

rh = m + E; , (4.16a) 
E = \ XX m ^4 - E m A'm B /r AB ) , (4.16b) 

1 A B 

1 2 A-\j2 m B/r AB ) =0, (4.16c) 

(4.16d) 
(4.16e) 
(4.16f) 
(4.16g) 

where to = J2A m A, and Q = Q n . In Eq. (|4.15|) we show schematically the e order and the 
r dependence of the terms neglected. Note that, by virtue of the Newtonian equations of 
motion, E and J 1 are constant to leading order. Here to, and J 1 are to be evaluated at 
u = t — R. Combining Eqs. (|4.15|) , fl3.5|), (|4.14|) and ( |4.13| ), we find, to the required order 
that P s j = -{7/6)mKQ ij . 

Combining Iq, Ip and Ig , we obtain finally 

j ew = Y. m ^i f 1 + \ v a ~ o Er^l + \Y. m AX A v\ 

A \ 2 2 B r AB J B A 

+4 l 2 ^ - H4 - 22v A • v B 

12 ^ r AB I 

-(v A • n^i?) 2 + 2(v B • n AB ) 2 - 2\ A ■ n AB v B ■ h AB 

x TO(7 1 

+2(a A + a B ) • x^b + 6 > f 

_ 1 1 1 [( ^ + _ (( ^ + Vfl) fi ^ )2]x ^ 

t^ AS ^ ^ 

+2(v A + v fl ) • x AB (10vM + IIuVb) - (26^ - ^M)r 2 AB 

~TK^l m Am B r AB \& A ■ h AB x { - A K rv'l B - a { }x A + 23 (a A + a B ) {l x A \ 
12 AS 

14 

-3^mM B + - —m-RjQ* + 0(e 3 ) x Q« , (4.17) 



AB 



where Q 1 ^ is a complicated 3-body term arising from the Pk m U,km term in Eq. (|4.8|), that 
vanishes identically for two-body systems. It is evaluated in Appendix D. 
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C. The three-index moment 



Since T$ w is dominantly of 1/2PN order, we need to calculate only the first post- 
Newtonian corrections to it, i.e. terms of 3/2PN order. We first note that P^ w [Eqs. 
( gl9bD and ( gjOg )] can be written 



jijk 
1 EW 



jijk , jjik fkij 
1 EW + 1 EW ~ 1 



EW > 



(4.18) 



where we separate I'ew into compact, field and surface contributions, given by 



7-iifc . fiife 



r 0i x j x k d 3 x , 



(d/dt)/ 



(1/16tt) / A H h jkl K. A d 2 tt . 

JdM 



(4.19a) 



Substituting Eqs. Q, ([pa]), (|3~8dl) , (gj|) and (gj) into (-g)T 0i and expanding through 
0(pe 3//2 ), we obtain 



jg* = £ m^s? ( 1 + k + 3 £ — J + O(e^) x 0* . 

A \ 2 B T ABj 



(4.20) 



To the required order, 



A 



16 



U, k (U kti -U itk ) + ~UUi 



(4.21) 



We then calculate Pp 1 following the method laid out in Sec. [IV A| . In the course of this 
calculation we find no TT terms dependent on positive powers of TZ. Finally we evaluate 
the surface contribution using Eqs. ( [4.14|) and (|4.15|) evaluated to lowest order, and find no 
contributions. The final result is 



1 ew = £ m * v a x a U + 

A \ B 

1 x m A m B „ ^ jk 
-^L — W A ■ n AB n AB x J A 



m B 
tab 



ab 



tab 



AB 



— m A m B r AB 2v A • h AB fil B + ll(2^n^ B - v A h% - v k A n 3 AB 



AB 



+^m A m B VA -h AB n%x h 2 + 7v A x ( in k 2 B -7v]ix k 2n t AB + 0(e 5/2 ) x Q ij . (4.22) 



The result agrees with Eq. (A52) of [ 35[ . The full moment can be constructed from 
this using Eq. ( |4.18|) . 



D. The four-index moment 



Since this moment contributes to the waveform already at PN order, we only need to 
evaluate the integrands through their first PN corrections. We write Eq. ( |2.19c|) in the form 
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jijkl _ jijkl , jijkl / » OQ \ 

(there is no surface contribution). Expanding {—g)T 13 through 0(pe 2 ), we obtain 

If 1 = £ m A vtx k i f 1 + + 3 E — ) + 0(6 3 ) x Q« . (4.24) 
To the required order, A tj can be written in terms of the instantaneous potentials 



A ij = AiUiUj - ^SijUkUk) + 16 



1 u ,(i X J) + 2U ,d u j) ~ u k,iU k j - U itk U jtk + 2U k;{i U j)>k 



1 •• 3 • 

-$ij(gU ik X jk + -U + U :k U k — U mik U[ m:k y) 



+ 0(pe 3 ) . (4.25) 



The term proportional to S tJ produces no TT contributions to the waveform, so we drop it. 

The method proceeds as in the previous cases, without the complications of cubic non- 
linearities. The result is 



Fiw = E m A (vl ~ \ E ^iBni B /r AB ] 4' + ^E rn A m B r AB n i i B (h'Z B - 5 kl ) 
1 



12 

A \ B / AB 



+ ^Y. m A<V i ix A l 
A A 



7 E m A m B% k A l r AB (?v l i + (Av 2 AB - v 2 B )n{ B 



4 

^ AB 



-3(v A • YL AB fn 3 AB - 4(3v A • n Ai? - 4v B • h AB )v { A n 3 l B 
— lQa B x 3 J B — 2a A x 3 J B + a A ■ ^i- AB n AB 



l AB 



-2j2 m c(l/rAc + l/r BC )h l j 
c / 

~h E ({Av AB - v 2 A )n% + 2(Sv i i B - 23vg) 

-(v A • n AB ) 2 h l { B - 8v AB ■ h AB v { AB n J 2 B + 4v A ■ n A B^AB 

-\Aa { Xx 3 l B + a A ■ x AB nl B - 4^(m c •, h 'ac)^ab 

c I 

E m Am B r AB ni B Uv 2 AB - 5v 2 A + 9(v A • h AB ) 2 - 4 E «W?"ac - 3a A • x AB ) 

+ - E m A m Bnis :c A ( V A^AB + Wv A + 4 (v A • n.AB ~ 3v A B • nAB^A^AB 
6 AB 

-3(v A • h AB ) 2 n l l B - \Aa ( A x 3 \ B + a A • x AB n l j B ) 

+ T^^2 m A , m B r AB nj LB (y A + 2v A • fi AB v A 7i AB - a { A x\ B + 2ajx5) 
iZ AB 

+ TK^l m Am B r AB n k AB Uv % i B - 2lv % i + 8v A b • nAB^AB^AB 

[Z AB 
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-6v A • nABVAn J } B + 35a^a^ B 

J AB n AB v AB n AB v A n AB v A n AB 



+ -^2m A m B r AB [2v 



3 as 



+ o Z m ^B \2v\n%v\ x A -v A - n AB ni B v\ x> A - l2v\n A ' B v\ B x A 



3 All 



■—mllQ ij 8 kl + 0(e 3 ) x . (4.26) 
oO 



E. The five- and six-index moments I B ^ m and I l B y/ nn 

These moments contribute to the waveform at 3/2PN and 2PN order, respectively, thus 
we only need to evaluate the dominant, Newtonian contributions to the integrands. Splitting 
the moments into a compact and a field piece, substituting the lowest order contributions 
to r ij at O(pe), into Eq. ( |T9g) , namely (-g)T ij = Ea™a^5 3 (x - x A ), and A ij = 



^(UjUj — \8ijU t kU t k), and carrying out the integration procedures as above, we obtain 

jijkim _ 1 d J ^ klm I ij 1 x - m B i:j \ 
Iew ~3lt\^ mAXA [ A 2^^ nAB ) 

4EwiA4(4-^)l+0(^) x Q« , (4.27a) 

72 



jijklmn _ ■'• ^ / fc/mn | ij ^ \~* r ^B^ij 

Iew ~vidTA 1 t mAXA Y A ~~^7 A ~ B nAB , 



X 

10 



+^ Z m Am B r AB nl B x { A (h AB > - 5 mn) ) 

1 AB L 

-x% {2h k A l ^ n - 2h { AB 5 mn) - 5 {kl 5 mn) ) 

n Qi3 6 {U 5 mn)^ + x gij _ (4.27b) 
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Equation (|4 . 2 7a|) is equivalent to Eq. (A53d) of p5 



V. EVALUATION OF RADIATION-ZONE CONTRIBUTIONS 

We now turn to the evaluation of the contribution h l ^_j^(t, x) given by the integral over 
the remainder of the past light cone of the observer, Eq. (|2.23|) . There is no material 
source now, so r u = A u /167r. On the other hand, the time dependence in the integrand 
of Eq. ( |2.23|) is not the simple fixed retarded time u = t — R of the EW moments. The 
(u' + r')-dependence of r u in Eq. Q2.23| ) reflects the variation in retarded time along each 
two-dimensional intersection of the past light cone of the event (t, x) with the future light 
cone of the event (w',0). However, r u is a functional of retarded potentials, such as V. 
When evaluated at u' + r', V has the form 
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V(u' + r',x')= f - f X " . <T(u' + r'-\x!-x"\,x."), (5.1) 
J \x — x'\ 

Notice that, because |x"| <C TZ, while |x'| > 1Z, we can approximate 

u' + r' - |x' - x"| « «' + n' ■ x" + 

+ (2r')" 1 [(n . x") 2 - r" 2 ] + . . . , (5.2) 

where n' = x'/r', and then expand such retarded functions about w' in powers of the small 
quantity r" / r'. For a given w', the retarded fields that contribute to A 1 - 7 along the intersection 
between the two light cones in Fig. 5 all have their source in the near zone, on slices of the 
near zone world tube that pass through the center of mass at time u' . The expansion ( |5.2| ) 
simply reflects the fact that, as one moves around the source in angle (integration over d 2 Q 
in Eq. ( |2.23| )), the orientation of the slice of the near-zone world tube that generates the 
fields precesses (see Fig. 6). 

Since the ingredients of A 1 - 7 are all fields evaluated in the radiation zone, we can use 
expansions in powers of 1/r', such as those of Eq. fl2.14p . The angular dependence of 
such expansions can always be expressed in terms of STF products of radial unit vectors n' 
(analogues of spherical harmonics). Thus A 1 - 7 can be written, in the regime r' ^> 71, as a 
sequence of terms of the generic form f l ^ l {u')n' <L> r'~ N . Then a change of variables 

C = {t - u')/r = 1 + (u - u')/r (5.3) 

puts Eq. ( |2.26| ) into the form, for each (N, I) term, 



/2\ N ' 2 r 1 + 2n / r dC 
W-m(N, I) = (-) I —-±— M u - r( C - 1)) 



x 



(C 2 - 1) 

(An)- 1 P d(j)' f 1 n' <L> (C ~ n ■ n^dcostf' 



N-2 



+ ';J W(c^TF^ fe(l '- r(c - 1)) 



x(4tt)- 1 d<p' ^ h' <L> ((- n' -h) N - 3 d cos 6' , (5.4) 



where (cf. Eq. ( |2.25| )) a = (( — l)(( + l — 2TZ/r)(r/2Tl). We first carry out the angular inte- 
grals, which yield n <L> An^((,-, a), where A^j can be computed from Legendre polynomials 
Pi(z) by ^(Ca) = \Sl- a Pi(z)(C ~ z) N ~*dz [see Appendix A, Eq. (^5j); a = 2 corre- 
sponds to the full 4ir angular integration]. Then, in the ^-integration from 1 to 1 + 2TZ/r, 
we expand the retarded time dependence of the f^j about u, then integrate; this is valid 
since 1Z < r. In the integrals from 1 + 21Z/r to oo, we integrate by parts numerous times, 
each time increasing the number of time derivatives of Jn,i, stopping when the result ex- 
ceeds the PN order required. The boundary terms that arise are evaluated at £ = 1 + 271 jr 
and ( = oo, corresponding to retarded time u — 21Z and — oo respectively. At the former 
boundary, we again expand the functions about u; at the latter boundary the contributions 
are assumed to be zero, which is equivalent to making the usual and reasonable assumption 
that the source is not extraordinarily dynamical in the infinite past. 
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For the cases where the field point is inside the near-zone, Eq. (|5.4j) still applies, except 
that now r < R, and the first ( integral runs from — 1 + 271/ r to 1 + 21Z/r (Fig. 5). 

In working to 2PN order, just as in the case of the EW surface integrals, Eqs. ( [2.20 ) 



here we must also evaluate the integrand A 4 - 5 to 0(pe 3 ). Here, as before, it can be shown 
that only the twice iterated fields are needed in practice. This can be seen as follows. We 
are only interested in the 1/r part of the waveform. Thus a contribution to A 1 - 7 ' that is 
already 0(pe 3 ) but that falls off faster than r'~ 3 (N > 3) can be dropped. This would apply 
to all terms that are quartically non-linear and higher, such as terms of the form h 2 (h t \) 2 , 
which fall off as r' -6 . Cubically non-linear terms of the form h(h y \) 2 can also be dropped; at 
leading order, they are 0(pe 2 ), but fall off as r'~ 5 . One might worry that by expanding 
in Eq. ( |5.4|) about u (the value of retarded time at which all contributions to the waveform 
are to be evaluated in the end), one could reduce the rate of fall off by one power of r 
for each retarded time-derivative. But each time-derivative either raises the order of the 
term by e 1 / 2 or kills it outright via a conservation law, such as for the Newtonian potential 
h ~ m/r. Thus the leading cubically non- linear contribution turns out to be of order pe 3 /r' 5 , 
which can be dropped. Thus only quadratically non-linear terms of the form (h y \) 2 need to 
be considered. As before, a knowledge of h a/3 to the accuracy shown in Eq. ( |3.5|) suffices; 
higher-order terms contribute terms of order pe 3 /r' 4 . However, we must now be careful in 
evaluating the terms which do contribute. For example a term of O(pe) that falls off as 
r' -6 , can, after three terms in the Taylor expansion of its retarded time dependence about 
u in powers of r(£ — 1), lead to a 1/r contribution to the waveform at 0(D _1 pe 5//2 ), which 
is 3/2PN order beyond quadrupole order. A term of this form would arise from the cross 
term between the gradient of the Newtonian potential m/r and the Newtonian quadrupole 
potential ~ /r 3 . Similarly a (pe)r' -7 term would contribute a 2PN contribution to the 
waveform. Such a term would arise from a cross term between the Newtonian potential and 
the Newtonian octupole potential ~ Q^ k /r 4 . A consequence of these considerations is that, 
in expanding the second- iterated fields h al3 , we must use the general multipole expansions 
of Eq. ( |2.14| ), expanded through octupole order. This amounts to expanding h 00 through 
q = 3, h 0t through q = 2, and h 00 through q — 1. Evaluating the integrals M al3kl '" kq to 
the needed order, using the general method for integrating over the near-zone hypersurface 
Ai described in Sec. |TD|, and adding any contributions to h a ^ from the radiation-zone 



integrations (primarily from W %3 ; see Appendix C), we obtain 

h 00 = Am/r' + 7(m/r') 2 + 2 tr'~ l Q ij (u')\ - - {r'~ l Q ijk {v!)) , 

l J ,ij 3 l J ,ijk 

h 0i = -2 {r'^lQ^iu') - e i3a J a (u')}\ + \ \r'~ X \Q iik {v!) - 2e ika J aj (u')]\ , 

h? = (m/r') 2 h' i3 + 2Q ij (u')/r' - - {r'~\Q ijk (u') - 4e (i|fca j a|j) («')]} fc , (5.5a) 



where Q l \ J a and J aj are defined in Eqs. fl4.16d| ) - ( 4.16g ), and where the superscript 



notation w a - k w denotes symmetrization only on i and j. 
To the required order, we then have 

A « = _ h wft3 + hi^X , + \h m ^h kk d) + 2h 00 H»° , (5.6) 

with the result 
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The terms in Eq. ( [5.7D are of the generic form f n ,iW)h' <L> r'~ N . We substitute each such 
term into Eq. Q5.4|), integrate using the procedure outlined above, and keep only terms 
through 2PN order that fall-off as 1/r. Evaluating at the detector distance R, we obtain, 
finally 
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As with the calculation of EW moments, we discard terms that fall off with increasing 1Z. 

The integrals involving the logarithm of retarded time are the tail terms, and are in 
complete agreement with [|3"3[1 , including the constants (11/12, 97/60, 7/6) added to the 
logarithms. Their origin is the backscatter of the outgoing gravitational waves off the lowest- 
order, Schwarzschild-like, static background curvature of the spacetime surrounding the 
source. More precisely, the logarithmic integrals can be seen to arise directly from the term 
—hP°h^ in Eq. (|5.6|) , which represents a modification of the flat spacetime characteristics 
by the potential h 00 ~ m/r. The first tail term, arising from the 2d A Q < ^ > j 'du 4 term in 
Eq. ( |5.7| ), is actually of 3/2PN order, while the second and third terms, arising from the 
ld 5 Q^ k /du 5 and ld 4 J aj /du 4 terms in Eq. (]577|), are of 2PN order. On the other hand, 
only the 3/2PN tail term contributes to the energy flux at 2PN order, resulting in the "47r" 
term for circular orbits in Eq. (|1.4j) . Notice that the tail terms show no dependence on the 
near-zone boundary radius 1Z. In the BDI framework, the tail terms contain a scale b which 
is associated with a gauge transformation from harmonic coordinates to a set of radiative 
coordinates used in that framework; physical results in the end do not depend on b, and the 
tail effects in the two frameworks are in complete agreement. 

It is easy to see that the final term in Eq. ( |5.8|) , which depends linearly on TZ exactly 
cancels the sum of the corresponding terms arising from the two- four- and six-index EW 
moments [Eqs. (pTp , flO§ and (jjTJj) ]. 

Thus combining the contributions of the six EW moments to Eq. ( |2.18| ) with these 
contributions gives the gravitational waveform, valid to 2PN order, for a general N-body 
system. The waveform is explicitly finite, with no divergent integrals or undefined terms. 
Henceforth, we shall not display any ^-dependent terms. 



VI. REDUCTION TO TWO-BODY SYSTEMS 

A. Center of mass and equations of motion to 2PN order 

We now specialize to the case of two bodies. Through 2PN order the dynamics of 
two-body systems are well known. The motion is governed by a Lagrangian that admits 
a conserved total energy and angular momentum, as well as a "conserved" center-of-mass 
definition. We define the system's center of mass X and the relative position x by 

X = m-^miX! + m 2 x 2 ) + f«(x 1 ,x 2 ) + f( 2 )( Xl ,x 2 ) + 0(e 3 ) x X , (6.1a) 
x = xi - x 2 , (6.1b) 

where m = m 1 + m 2 , and and denote 1PN and 2PN corrections to the center-of-mass 
definition. Inverting these expressions and setting X = 0, we obtain 

Xl = (m 2 /m)x - f (1) - f (2) + 0(e 3 ) x Xl , (6.2a) 
x 2 = -(mi/m)x - f (1) - f (2) + 0(e 3 ) x x 2 . (6.2b) 
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The only place the 2PN correction f ( 2 ) could conceivably be needed is in the lowest-order 
quadrupole moment, but in this case it is straightforward to show that it is not, since 



Q ij = Y, m AX l A x j A = /ixV + m/ (1) V (1)i + 0(e 3 ) x Q ij , 



(6.3) 



where f^ 1 = ~^i](5m/m)(v 2 — m/r)x l (see, e.g. fl59|), and where we define the two-body 
variables fi = 17111712/171 (reduced mass), 77 = fi/m, 5m = m\ — m2, v = Vi — V2, and r = |x|. 
The two-body equations of motion then take the effective one-body relative form, through 
2PN order: 



a — a^v + a plr + a 5o ^ a 2pV "I" a ss 0(&^^ ^) 



(6.4) 



where the subscripts denote the nature of the term, post-Newtonian (PN), spin-orbit (SO), 
post-post-Newtonian (2PN), and spin-spin (SS); and the superscripts denote the order in e. 
The individual terms (excluding spins) are given by 
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B. Two-body Epstein- Wagoner Moments 

Restricting the summations in the EW moments to two bodies and substituting Eqs. 
>72|) , we obtain, through 2PN order, 
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In addition, the moments that appear in the tail terms, Eq. 
order, to 



(6.6c) 
(6.6d) 
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|), reduce, to the required 



Q ij = fix ij , 

Qv k = -^(5m/m)x l3k , 
J aj = -fi(Sm/m)(x x \) a x j . 



(6.7) 
(6.8) 
(6.9) 



C. Two-body gravitational waveform and energy loss 



Substituting the EW 2-body moments ( |6.6|) into Eq. (|2.18|) , calculating the time deriva- 
tives using the 2PN equations of motion ( |6.5| ) to the accuracy needed, and adding the tail 
terms from the radiation zone integral (|5.8| ), we obtain the gravitational waveform. An 
alternative method is first to calculate the so-called "symmetric trace-free" (STF) moments 



defined by Thorne |34[ and used by BDI, and then to calculate the waveform. The proce- 
dures and formulae needed to do this are given in Appendix E. The result for the waveform 
is 
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where, as before, the superscripts denote the effective PN order, and subscripts label the 
nature of the term, and where the individual non-spin pieces are given by 
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The leading PN and 3/2PN spin-orbit and the 2PN spin-spin contributions to the waveform 



can be found in Eqs. (3.22) of [41| and in Appendix F. There will also be in principle 2PN 
spin-orbit terms; these have not been calculated to date. 

Although we have differentiated the moments appearing in the tail terms explicitly using 
the equations of motion in order to display the waveform contributions in a consistent 
manner, this is not the best form of the tail terms for explicit numerical evaluation in the 
case of general orbits. The reason is the slow fall-off of the logarithmic term with increasing 
s. Instead, it is preferable to revert to the forms of the tail terms given in Eq. ( |5.8| ), split 
each integral over s into a finite part from to sq, where so corresponds to several dynamical 
timescales of the source, and a remaining integral from so to oo. The first integral can be 
done using the expressions given in Eqs. (|6.11|) . The remaining integral is integrated by 
parts twice. One can then show |J3] that the latter integral falls off as 1/ So generally, and for 
nearly periodic orbits, as 1/s 2 ,. By choosing sq sufficiently large (generally a few dynamical 
timescales or orbital periods), one then can obtain accurate numerical representations of the 
tail terms, without having to integrate over the entire past history of the source. 

Differentiating h lJ with respect to time, using the 2PN equation of motion (|6.5| ) where 
required, and substituting into Eq. ( |2.30|) ; or equivalently, taking the appropriate time 
derivatives of the STF moments (Appendix E), and substituting into Eq. (E5b| ) , one finds 
for the energy flux, 
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where the non-spin contributions are 
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The 3/2PN spin-orbit and 2PN spin-spin contributions can be found in Eqs. (3.25) of [41 



and Appendix F. The tail contribution is formally of 3/2PN order, arising from a cross term 
involving P 3 ^ 2 Qtail an d U ! f° r simplicity, we do not write it out explicitly (for circular 
orbits we evaluate it below). The "11/12" term in Eq. Q5.8| ) contributes a term of the 
schematic form (d A Q/du A )(d 3 Q/du 3 ), which can be written as a total time derivative and 
absorbed into a redefinition of the energy E at an order above that at which it is well 
defined as a conserved quantity (see e.g. |60|j61[1 for a discussion of this point). In the same 
way, the form of the tail term shown in Eq. (|6.13cj ) has been achieved by integrating the 
tail contribution once by parts and moving the total time derivative over to the left-hand 
side. The 2PN tail terms in the waveform make no contribution to the energy flux to 2PN 
order because their cross product with the quadrupole piece contains an odd number of unit 
vectors N, and thus vanishes on integration over solid angle. They will, however, produce 
5/2PN contributions to E via cross terms with the 1/2PN waveform terms P x l 2 Q l K 
Through first PN order, Eqs. (|6.13|) agree with |p!7| , |62| . 



VII. QUASI-CIRCULAR ORBITS 



A. Orbit equations and gravitational waveforms 



Because gravitational radiation reaction circularizes orbits, the late stage of inspiral of 
a compact binary, such as that of the binary pulsar PSR 1913+16, will be characterized by 
a quasi-circular orbit, that is, an orbit which is circular apart from the slow inspiral caused 
by radiation damping. We define the Newtonian angular momentum L^r = /jx x v, the unit 
vector A = Ln x n, and the angular velocity uj = |LAr|//jr 2 . A circular orbit is given by 
the conditions r = r = 0. Solving the 2PN two-body equations of motion ( p. 5| ) under these 
conditions gives 



UJ' 
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Then the orbital velocity is v = ruX and the orbital energy through 2PN order is 

\m \ 1 /77i x 2 
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In order to calculate waveforms as observed by an Earth-bound detector, we must choose 
conventions for the direction and orientation of the orbit. The standard convention is to 
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choose a triad of vectors composed of N, the radial direction to the observer, p, lying along 
the intersection of the orbital plane with the plane of the sky (line of nodes), and q = N x p 
(see Fig. 7). The normal to the orbit L^r is inclined an angle i relative to N (0 < i < tt). 
The orbital phase — ^ u + const of body 1 is measured from the line of nodes in a positive 
(out of the plane) sense (orbits seen to be moving clockwise correspond to i > vr/2). The 
two basic waveform polarizations h + and h x are given by 

h+ = ^(PiPj - QiQjW 3 , (7.3a) 
hx = ^(PiQj + qiPj)h l] . (7.3b) 

(There is no need to apply the TT projection in Eq. (|6.10| ) before contracting on p 
and q.) From our conventions, we have that n = pcos0 + (qcosi + Nsini)sin0 and 
A = — psin0 + (qcosi + Nsinz) cos0. Since h % i consists of terms of the form h % h\ A*A J or 
n^X^, we find the following formulae to be useful in evaluating the polarizations: 

(n i n j ) + = - sm 2 i + -(1 + cos 2 i) cos 20 , (7.4a) 
1 1 

(A l A J )+ = -sin 2 i- -(1 + cos 2 i) cos20 , (7.4b) 

(n( i \^) + = -1(1 + cos 2 i) sin 20 , (7.4c) 

(n l h^) x = - cos i sin 20 , (7-4d) 

(A i A i ) x = --= cos i sin 20 , (7.4e) 

(h {i \ j) )y, = ^ cos i cos 20 , (7.4f) 

N ■ n = sin i sin , (7.4g) 
N-A = sini cos0 . (7.4h) 

Substituting r = and Eq. ( |7. 1|) into Eqs. (|6.11| ) (keeping PN and 2PN corrections in Eq. 
( [i 7 . 1| ) as needed), and using Eqs. ([r.4[), we can evaluate h + and h x explicitly as functions of 
orbital phase and orbital orientation. The waveforms can be expressed in terms of powers 
of m/r, but it is observationally more useful to express them in terms of muj m (m/r) 3 / 2 , 
since u is directly related to the observed gravitational-wave frequency. Instead of showing 
the result here, we refer the reader to where the complete, "ready-to-use" pair of 2PN 



waveform polarizations are displayed and discussed. Similar substitution into Eqs. (|6.13|) 
results in Eq. ([T 



B. Tail Terms 

Because they involve integration over the past history of the source, the tail contributions 
to the waveform and energy flux require additional discussion. For circular orbits, the + 
and x polarizations of the quantity P 3 ^ 2 Qtail are given by 



36 



(P 3/2 Qtail)+ = 8m(l + cos 2 i) J ( — cos 2. , 



v 
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(7.5b) 



Because r and u evolve on a radiation-reaction timescale trr which is long compared to 
an orbital period, we can approximate them to be constant in the above integrals; the 



results will be valid up to corrections of order {ojtrr) -1 ^(ujtrr) <C 1 fl45f . Notice that the 
integrals converge as s — > oo, even if we approximate m 2 /r 4 « constant (in fact, r — > oo in 
the infinite past [|63|| , so the integrals truly converge). Thus we can substitute cj(it — s) for 
with u; = const in the tail integrals, pull out the m 2 /r 4 factor, and use the fact that, for 
any integer n, 



1 s = 



' c = 



sin(na;s) In 



cos(nus) In 



— ^ — ) ds = — -(7 + \n(2nuR) + 0[(2nuR)' 2 )) 



where 7 is Euler's constant. The result is 
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It is useful to combine these tail terms with the lowest-order quadrupole terms, given from 
Eq. ( S.lla|) by Q + = — (m/r)(l + cos 2 i)cos20 and Q x = — 2(m/r) cosi sin20, into the 
forms 
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where 



if) = - 2(m/r) 3 / 2 [ 7 + ln^i^- 11 / 12 )] 
= cu{u - 2m In/? - 2m[7 + ln(4cje" 11/12 )]} 



(7.9) 



We first note that one effect of the tail term is to shift the phase of the quadrupole piece 
by an irrelevant constant, and by a term which varies logarithmically with u as the inspiral 
proceeds. This slowly varying phase shift was studied in |38|| . 

We also recognize that u — 2m In R = t — R — 2m In R is retarded time with respect 
to the "true" null cone that intersects the observation point at (t,R). This can be seen 
by noting that, in the asymptotic, Schwarzschild-like spacetime of the source, in harmonic 
coordinates, outgoing radial null geodesies obey t — r — 2mlnr + 0(1/ r) = const. An 
identical .R-dependence in the phase shows up at the next 1/2PN order, when one combines 
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the two polarization states of P l / 2 Q^ with those of P 2 Qtail- We thus conclude that, 
at least through the 2PN order considered, our procedure for calculating the tail terms 
yields gravitational waves that asymptotically propagate along the true harmonic null cones, 
toward true future null infinity, despite the use of a flat-spacetime wave equation for h a/3 . 
This avoids the need for further matching or other devices to connect our solutions to true 
null infinity, and answers another long-standing criticism of the EW framework 0. It is 
useful to note also that, in the BDI approach, a similar logarithmic term appears in the 
phase shift ( [7.9|) , but there the term depends on the parameter b used in the transformation 
from harmonic to radiative coordinates. The appearance of such a parameter can be shown 
to have no physical consequences, as expected |3£|,|4j]. Our method is explicitly free of such 



arbitrary parameters, all effects of 1Z having cancelled. The only external radius which 
appears is that of the observer. 

The tail contribution to the energy flux, given by Eq. ( |6.13c|) can also be calculated in 
closed form using the above assumptions together with Eq. ( |7.6b| ). The result is the "47r" 
term in Eq. (|1.4|). 



C. Display of the waveforms 

We now display our results explicitly by plotting the waveform for an inspiralling binary 
as a function of time. We will assume that the binary is in a quasi-circular orbit in its 
last few moments before the final plunge to coalescence. The time evolution of the orbital 
phase-velocity in this regime can be obtained by integrating the equation 



dio 



E 



dE/du 



(7.10) 



where E is given by Eqs. (|1.4| ) and dE/du can be obtained from Eqs. ( |7. 1| ) and ( |7.2| ). The 
orbital phase angle can, in turn, be obtained by integrating the orbital phase velocity. 
The results are 
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-3/8 



(7.11b) 



where T is a dimensionless time variable related to the coordinate retarded time u by T = 
rj(u/5m), and C and T c are constants of integration. The constant T c is the dimensionless 
retarded time at coalescence (the time at which the frequency in Eq. ( |7.11|) formally becomes 
infinite), and <p c is the orbital phase at coalescence. 

We can now use the orbital phase evolution along with Eqs. ([7.3|), (|7.4| ) and 1| ) to 
write h + and h x as explicit functions of time. We will not display the result here (there 
are enough large equations in this paper already), but rather refer the reader to Eqs. (2) to 
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(4) in |^8| for "ready-to-use" waveforms. The "ready-to-use" waveforms are essentially Eqs. 
( |6.11| ) boiled down to the circular orbit case. 

For the case of a 1.4M Q neutron star spiralling into a 1OM black hole the resulting 
frequency sweep and waveform are shown in Fig 8. The observer is viewing the orbital 
motion edge on, so that i = 71/ 2 in Eqs. ( |7.4| ). In this case the gravitational radiation is 
linearly polarized (only h + is present). The upper cut-off frequency in Fig. 8 is chosen to 
be 180 Hz; this is approximately the orbital frequency at the innermost stable circular orbit 
|65| , |66[| for this type of system. For the initial LIGO detector, Finn 67] has shown that a 



substantial fraction of the signal-to-noise ratio available is accumulated when integrating a 
matched filter against the signal in the frequency range we have displayed. In other words, 
the segment of the waveform shown in Fig. 8b, sweeping from 75 Hz to 180 Hz, is the 
portion of the waveform which is actually most detectable for the initial LIGO detector. 

As energy is extracted from the system by the radiation, the orbital radius shrinks and the 
orbital frequency increases. This gives rise to the dominant "chirp" feature of the waveform 
in Fig. 8b: the growing amplitude and the bunching of peaks at late time. However, 
because the coordinate velocity rises to about 0.5c, this system is quite relativistic, and thus 
the inclusion of higher multipoles of the radiation causes the waveform to differ considerably 
from the simple cosine chirp that one would compute just using quadrupole radiation. The 
pairing of wave crests (alternately closer together and farther apart) signifies the onset of 
the gravitational analogue of synchrotron spikes. Just as in electricity and magnetism this 
feature comes from the inclusion of many harmonics of the radiation. In our analysis we 
have included multipoles through the six-index multipole l l ^J nn . This allows us consistently 
to include components of the radiation in our waveform at multiples of the orbital phase 
n0 or bitai where n ranges from 1 to 6, n = 2 being the dominant quadrupole contribution. 

Another interesting feature of Fig. 8b is that adjacent troughs are not the same depth, 
but adjacent crests are essentially the same height. This effect also has a discernable physical 
origin. The deeper troughs arise when the lighter mass is coming toward the observer; thus 
the observer is in the forward synchrotron beam pattern of the lighter, faster-moving mass. 
The shallower troughs arise when the lighter mass is receding from the observer. [At the 
left-hand-side of the figure, the phase is arbitrarily set to zero, i.e. the heavier mass (chosen 
to be mi) is passing through the ascending node coming toward the observer and the lighter 
mass is receding. The waveform is clearly in the not-so-deep trough at this left-most point.] 
The crests are essentially the same height because the radiation is virtually the same when 
the masses are moving transverse to the line of sight of the observer regardless of which mass 



is closer to the observer (see [35] for further discussion of the asymmetric radiation emission) 



The extent to which the harmonic structure might be measurable by a gravitational-wave 
detector is currently under investigation |68| . Preliminary analysis shows that neglecting the 
harmonic structure {i.e. just using the quadrupole amplitude to describe the wave) results 
in approximately a 4% loss in signal-to-noise ratio. In Appendix F we show how the effects 
of spin modify the waveform and frequency evolution. 



VIII. DISCUSSION 

We have extended the Epstein- Wagoner framework for calculating gravitational radiation 
from slow-motion systems to produce a method that is free of divergences or undefined inte- 
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grals. The extension involved adding to the original framework the integral of the effective 
source over that part of the past null cone of the field point that is exterior to the near zone. 
When expressed in appropriate variables, that integral can be shown to be convergent, and 
can be evaluated in a straightforward way, to any chosen PN order. The exterior integral 
yielded (a) terms that explicitly cancel terms from the EW framework previously thought 
to be divergent (b) tail terms, in agreement with other methods based on matching, and 
(c) phasing terms that verify that the radiation asymptotically propagates along true null 
cones of the curved spacetime. 

This new, well-defined framework, provides a basis for extending the calculation of grav- 
itational radiation to higher PN orders. An extension to 5/2PN order in the BDI framework 



has been achieved by Blanchet ||69|| ; such an extension in the improved EW framework is 
in progress. Extension to 3PN order will be a bigger challenge, simply because of the com- 
plexity of the terms, including quadratically non-linear integrals, and the rapidly increasing 
number of computations. However, we foresee no obstacle in principle to such an extension 
in the improved framework. 

This improved framework will also allow derivation of near-zone gravitational fields in 
a form that will yield equations of motion for the sources to high PN orders. It should be 
possible to derive radiation-reaction terms in the two-body equations of motion, at order e 5//2 
and e 7 / 2 beyond Newtonian gravity [ 30, 31, 70 1, without the presence of ill-defined or divergent 
terms, and without the need for matching between zones. One goal would be to derive the 
non-dissipative, 3PN terms in the equations of motion. This would improve the accuracy 
of estimates, using a hybrid Schwarzschild-PN equation of motion, of the transition point 
between inspiral and unstable plunge in the late stage of compact binary inspiral p5| , |66 



Calculation of the near-zone fields will also be important in developing interfaces between 
the post-Newtonian approach, which works well for most of the inspiral, and numerical 
relativity methods which must be used for the final few orbits and the coalescence. Work 
on this latter subject is in progress. 
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APPENDIX A: STF TENSORS AND THEIR PROPERTIES 

In calculating field integrals we make frequent use of the properties of symmetric, trace- 
free (STF) products of unit vectors. The general formula for such STF products is 
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<L> 



p=0 



(21- 1)!! 



(Al) 
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where [1/2] denotes the integer just less than or equal to 1/2, the capitalized superscripts 
denote the dimensionality, I — 2p or p, of products of n l or 8^ respectively, and "sym(q)" 
denotes all distinct terms arising from permutations of indices, where q = l\/[(2 p p\(l — 2p)\] 
is the total number of such terms (see p4] , p8| for compendia of formulae). For convenience, 
we display the first several examples explicitly 
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There is a close connection between these STF tensors and spherical harmonics. For example, 
it is straightforward to show that, for any unit vector N, the contraction of N L with n <L> 
is given by 
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where Pi is a Legendre polynomial. This latter property can be used to establish the identity 
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Since the left- hand- side is STF, and depends only on the unit vector n, then it must be 
proportional to the STF combination n <L > . To establish the normalization, contract both 
sides with the //-dimensional non-STF product N L , where N represents the z-direction. 
Using Eq. ( |A3[ ), and recalling that Py = A/j'^/'o, where Mv is a normalization coefficient, 
we find that the integral yields [/!/(2Z — 1)!!]P^(N ■ n)5;/', establishing the unit coefficient in 
Eqs. ( PIP and (fA|). 

In calculating the radiation-zone contributions to h l \ we must also evaluate the integrals 
/ 2 " d(j)'jl_ a n' <L> (( - n' • h) N - 3 d cos 6', where a = (C~ 1)(C + 1 - 2K/r)(r/2K). The 
result must be an /-dimensional STF tensor, dependent on the only vector in the problem, 
n, and thus must be proportional to n <L> . To determine the proportionality factor, which 
will be a function of ( and a, we contract with n L , chose n to be in the z-direction. and 
substitute Eq. ( |A3D . The result is 
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APPENDIX B: DERIVATIVES OF GRAVITATIONAL POTENTIALS 



In evaluating the "field" parts of EW moments, we have repeated occasion to integrate 
expressions involving two derivatives, spatial, time, and mixed, of the potential U and two 
spatial derivatives of X. For a field point external to the bodies, such derivatives can be 
calculated easily from the expressions ( [4.4a] ) and ( |4.4b| ). However, because the integrations 
run over the locations of the bodies themselves, we must carefully evaluate the singular 
behavior of such double derivatives at x = x^. Consider, for example, the expression for U, 
written in terms of a smooth density distribution: 
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where a' = dv'/dt. For a field point outside the bodies, shrinking the density distribution to 
a point yields a result equivalent to that obtained by differentiating Eq. (|4.4a| ). For a point 
inside, say, body A, we find that the integral / body A Ud 3 x — > —(4'K/3)mAv\ as the size of 
body A shrinks to a point. Consequently we must add a delta- function term to all double 
derivatives of U and X found using Eqs. ( 4.4a ) and ( f4.4b| ). The results are 
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where f denotes derivatives computed from Eqs. (|4.4a ) and ( |4.4b| ) 



APPENDIX C: THE SECOND-ITERATED FIELDS 
In Section IIII BL we wrote down the second-iterated solution for h a/3 in terms of the 



potentials V, Vi and Wy. Here we discuss the solutions for these potentials, Eqs. (|3.4j ), in 
more detail, especially the potential Wy, whose source is non-compact. 

We first consider field points in the radiation zone. Since their sources have compact 
support, the potentials V and Vi do not have to be divided into contributions from integrals 
over the near zone and over the radiation zone. They can be expanded using the analogue 
of Eq. (|2.14j) , and written to the needed order in the form 



V(f,x) =m/r + \ {r-'QKu)} .. ~ \ K^V)} .., 

+ 1 -Q/r-{r- 1 F\u)} i + 0(e 3 ), (Cla) 
V l {t^) = -\{r- l [Q^{u)-^r{u)]} ] 

+ 1 \r- 1 [Q iSk (u) - 2e ika J a \u)]\ + 0(e 5/2 ) , (Clb) 
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where Q and F l = J2a m Ax\ (y\ — J2b m B l^f ab) respectively represent the difference 
between the monopole and dipole moments of the potential V, and the 1PN accurate, 
constant total mass m, Eq. ( |4.16a| ), and the vanishing center of mass X, Eq. ( [4.16c ). In 



constructing h 00 using Eq. ( 3.5a| ) these two terms are cancelled by terms from W = W n 

The potential W^ must first be divided into near-zone and radiation-zone contributions, 
Wij = (Wij)_sf + (Wij)c-Af- To the 0(e 2 ) needed for the use of Wy in source terms for 
higher iterations [see Eqs. (|3.5| )1, we can approximate the integrand in Eq. ( |3 . 4cj ) by <Jij + 
(47r)~ 1 (f/ ii f/ J - \5ijU ]k U ik ) = r^/4, with = Y,A m AV A v A 5 3 (x - x A ). Here denotes 
the first-iterated effective stress-energy. Then (W^at can be expanded using the analogue 
of Eq. (pT4p , with the result 

(Wij)M = £ (±M^- k A , (C2) 



where 



M i]kl - k "(u)= [ T;Uu,x)x kl - k «d 3 x . (C3) 

>M ( ' 



Using the expression above for in each of the moments in Eq. (|C3l ), and using the 
strategy for evaluating field integrals described in Sections IV.A and B, we find, to the 
needed accuracy, 

M ij = ~Q ij , (C4a) 

M ijk = _Qijk _ l e (i\kaja\j) ^ ^ C4b ^ 

6 3 

M m = — m 2 K(25 i{k 5 l)j - -5 lj 5 kl ) 
15 2 

+term independent of 1Z , (C4c) 

where we discard terms that fall off with increasing 7Z, but retain all other terms. Although 
we never actually need the contribution from the moment M^ kl , we show the 7£-dependent 
term to illustrate its ultimate cancellation. 

To evaluate (Wij)c-j\f, we use the fact that, to the required order, in the radiation zone, 
Tqn = (47r) _1 (m 2 /r /4 )(n /<y> — Using Eq. (|5.4[) , and remembering the factor of 4 

difference between h 1 ^ and Wij, we obtain 



{W^c-m = - ~^ <v> , (C5) 



where we again discard terms that fall off with TZ. It is easy to see that the 7^-dependent term 
in Eq. (|C5|) exactly cancels the corresponding term in (Wij)jy resulting from Eq. ( |C4c| ) and 
(p2|). Combining the contributions to W^ through octupole order, and substituting them 
along with Eqs. ( |CTD for V and Vi into Eqs. ( |3.5|) yields the second-iterated radiation-zone 
fields h a/3 , Eqs. ( |5.5| ). It is interesting to note that the (m 2 /4r 2 )n i - J term in (Wij)c-j\f is 
required in order that the far-zone field correctly approximate the Schwarzschild geometry 
in harmonic coordinates in the static limit, namely 
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h 00 = Am/r + 7(m/r) 2 , (C6a) 
h 0i = , (C6b) 
h ij = (m/r) 2 n ij , (C6c) 

[compare Eq. ( |5.5| )1. This contribution could not have been found using the EW approach 
without our new formulation of the radiation-zone integrals. 

We next consider field points within the near zone. Expanding the retardation about 
t = u with |x — x'| as the small parameter, we obtain Eqs. ( |3.6|) and (|3.7p . The compact 
contributions to U, X, Ui and Pjj can be evaluated directly; the non-compact part of Pjj is 
left unevaluated until it is incorporated into an EW moment (see Appendix D). It remains 
to evaluate the radiation-zone contribution (Wij)c-j\r with a near-zone field point. Using 
the form of tJL above, and using the near-zone field-point version of Eq. ( |5.4|) , we find only 
contributions proportional to m 2 r 2 /TZ 4 and m 2 /TZ 2 . Thus we can discard such terms. 



APPENDIX D: CUBIC NON-LINEARITIES IN I% w 

At 2PN order, the non-linear field source A 00 Eq. (|4.8|) contains terms that are cubically 
nonlinear, i.e. that depend on effective products of three gravitational potentials. The 
contribution of the final such term in Eq. (|4.8|) , proportional to UU tk U^ k , to the integral 
J M A 00 x' l x^d 3 x can be evaluated straightforwardly by integrating by parts. However, the 
two terms 2P jk U jk — P km U^ km are more difficult because Py itself [Eq. (|3.7d|) l is a potential, 



one of whose pieces is produced by a non-linear source. The contribution of the compact 
source can be handled easily by the methods of Sec. [IV A| . Here we focus on the non-linear 
piece. We define the non-linear potential 

Pij (u,x) = T^^iUiUj - ^5 tJ U k U k )(u,x') , (Dl) 

We then need to evaluate the integral 

(1/tt) / (2p, k U k - p km U km )x i x J d 3 x . (D2) 

J M 

We integrate the first term by parts, show that the surface terms fall off with TZ, and obtain 
8 J2a m Ap{~x.A)x % A ~ (4/ tt) / pU'^x^d 3 x. The first of these terms may be evaluated using the 
non-linear pieces of Eq. ( |4.4dj) . The second term may be written in the form 



\VU'\ 2 d 3 x r 



... Z'nJ 1 ,|V^|f fe. (D3) 

2ty z Jm a Jm\x. — x'| |x — ~x.a\ 

In the x-integration, we change variables to y = x - and integrate using the general 
method described in Sec. [IV A| . The result is the integral (1/2%) J2a m A Jm \^'U'\ 2 d 3 x' (x n] — 



x l j/)/\x.A — x'|, which can be easily evaluated by integrating by parts 
The second term in Eq. (|D2|) can be written 
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j M {y[ k u[ m --5 km \v'u'\ 2 ) d 3 x' 

T, m J f 3(3 i " ~ ^ km 5 3 (* - *a)) ^d 3 x (D4) 

Jm x — x' V x — xa 3 / 
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Again we do the re-integration by changing variables to y = x — x^, and using the method 
of Sec. fV A . Integration of the delta-function term is straightforward. The remaining 
x'-integration takes the form 



-/ (u[ k U[ m - l -5 km \VU'f)d*x> 

TT J M V 2 / 

1 5kmX ij 1 \ 

~3 |x' - + X fik(i$3)™ ~ ^^k(i^j)m) , (D5) 

where <& A = |x' — x^^/lS, ty A = |x' — x^| 3 /3, and X A = |x' — x^|. The first five terms 
in Eq. ( p5| ) can be evaluated simply by integrating by parts. The sixth term is equivalent 
to the cubically nonlinear term in A 00 proportional to UU^U^ [Eq. (|4.8|) ]. The final term 
proportional to TZ is straightforward. 

The seventh term requires extra work. Dropping contributions with no TT part, we 
find that the integral to be evaluated is 7r _1 j M U^UjXd 3 x. Defining XJ a and Xa to be the 
contribution to U and X from body A, respectively, we write 

/ U ti U d X<Px = W U A4 U Aij X A d 3 x + W U A .iUA,jX B d 3 x 

JM A JM A _£ B JM 

+ 2 E / UA.(iU Bj) X A d 3 x + E / UA,( t UB, 3 )X c d 3 x . (D6) 
a+b Jm a±b+c Jm 

The first term has no TT part, while the second two terms can be evaluated using the 
standard methods of Sec. |1V A] , and lead to the term —3J2AB m A m Bn A B i n Eq. ( (4.17| ). 
We define the third term to be G%\, change variables to u = x — xc, y = — xc, and 
z = x^ — xc, verify that no surface contributions at TZ are so generated, and show that Q can 
be written Q 1 ^ = T,abc rn A m B rn c V l y V{F(y, z), where F(y, z) = J M |u-y|~ 1 |u-z|~ 1 'ud 3 'U. 
The latter step involves ensuring that the piece of F that diverges with TZ contributes no TT 
part to Q, so that the integration can effectively be commuted with the y- and ^-derivatives. 
Note that F has units of (distance) 2 , is symmetric on y and z, is a function only of |y|, 
|z| and w = |y — z|, and has the property that V 2 F = —Any/w, V 2 F = —4nz/w. It is 
then straightforward to show that the function with these properties is given by -F(y, z) = 
— (27r/3)[(y + z)w — yz + (y 2 + z 2 — w 2 ) ln(y + z + w)], modulo terms that give no TT 
contribution to Q. Thus the solution for QlL in Eq. (|4.17| ) is 

E m A m B m c V A V :i B F(x A c, *bc) , (D7a) 

A^B^C 
2 

-- [{r AC + r BC )r AB - tac^bc + 2x AC • x BC ln(r AC + r BC + r AB )\ . (D7b) 

Note that, because V l A V B {x A c • *bc) — no logarithmic dependence on source variables 
actually survives in hj, T . For two-body systems, this term does not enter the formula for 
the EW moment. 



G lJ — 
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APPENDIX E: STF-MULTIPOLE DECOMPOSITION 



Although the Epstein- Wagoner multipoles arose very naturally in our retarded-time ex- 
pansion of the relaxed Einstein equation, these are not the only multipoles for displaying 
the answer. An alternative set are the symmetric tracefree (STF) multipoles, which arise 
naturally in angular decompositions of the waveform (see, e.g. [0), and are multipoles of 
choice in the BDI framework. Thus it is useful to obtain a transformation between the EW 
multipoles and the STF multipoles. 

If the waveform is known then the STF multipoles can be projected out. This is exactly 
analogous to projecting the coefficients of spherical harmonics from a scalar function. The 
STF multipoles can be projected from the TT waveform by integrating over the sphere [see 
0, Eq. (4.11)]: 
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where Xg T l F "' am are called "mass" multipole moments and Jp r a p" am are called "current" or 
"spin" multipole moments. Substituting the expansion of h T ^ 2 in terms of EW moments, Eq. 
( |2.18| ), and adding the radiation- zone tail terms, Eq. ( |5.8| ), we obtain the transformations, 
correct to 2PN order: 
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where the STF notation on the right-hand side means symmetrize and remove all traces 
(note that the STF tensors are symmetric on all indices, while the EW moments are sym- 
metric only on selected pairs). These transformations can also be established using Eqs. 
(5.23) and (5.24) of @. 
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For two-body systems in general orbits, the resulting STF moments are given by 
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where the "tail" STF moments are given by 
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Through 3/2PN order, these moments agree with p5|, and in the circular orbit limit, 



through 2PN order, they agree with BDI [3£ 



In terms of STF moments, the waveform and energy flux may be written []34 
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Substitution of Eqs. (|E^) into Eqs. ( |E5a|) and (|E5b| ), using 2PN equations of motions in 
any acceleration terms generated by time derivatives, and keeping terms through 2PN order, 
yields Eqs. flPOD , ( |67TTD , ( |67T2l) and (|6TT3| ). 



APPENDIX F: SPIN EFFECTS 

In this paper, we have used our augmented Epstein- Wagoner formalism to give a complete 
description of the gravitational radiation for inspiralling "point-mass" binaries through 0(e 2 ) 
beyond the lowest-order quadrupole contribution. In this Appendix we demonstrate that our 
formalism is also adequate for computing contributions to the radiation which arise from the 
finite spatial extent of the bodies. Our primary goal will be to compute the contributions to 
the radiation from the bodies' spin angular-momenta, but in the process we will show how 
other extended-body effects, such as those due to a body's intrinsic quadrupole moment, 
could be computed with our formalism. The results will be presented in such a way that the 
spin contributions computed here can just be added to results already presented here and 
elsewhere. In particular we give the spin-orbit (PQ^o an d P^^ 2 Q % so) an d spin-spin (P 2 Q l ss) 
contributions to the waveform Eq. ( |6.10| ) for general orbits. We also give a restricted 



circular-orbit version of the results which can be added to the "ready-to-use" waveforms in 



In order to derive the spin corrections to the waveform, we relax our "point-mass" as- 
sumption and allow the bodies to have spatial extent small compared to the interbody 
distances. We further assume that the bodies are uniformly spinning fluid balls, approxi- 
mately spherical in harmonic coordinates. (A full discussion of this "fluid sphere" formalism 
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is given in Appendix A of ||35|| , where it is used to derive the waveform produced by non- 
spinning bodies through Ofe 3 / 2 ].) Although formally, our PN approach restricts us to weak 
internal gravity, we anticipate applying the results to neutron stars and black holes, as in 
the non-spinning case, by relying upon the Strong Equivalence Principle (see Sec. II. B for 
discussion of this point). It is now conventional, in treating spinning compact bodies, to 
view the spin S of each body as a quantity measured in units of its (mass) 2 , as is the case 
for black holes. Given that, formally, 5* ~ mdv, where d is the size of the body, and v is 
its rotational velocity, our convention implies that Scompact /^Formal ~ m/dv ~ e 1//2 , with the 
result that spin effects are viewed as 1/2PN order smaller per factor of spin than would be 



the case formally (see p^fflf for further discussion). 

The leading-order spin corrections to the waveform arise solely from terms in the source 
[Eq. (|2.5|)1 directly dependent upon fluid velocities. Since these terms have compact support, 
they generate no contributions to the waveform from surface terms or from far-zone integrals, 
at the order we are considering in this Appendix. Thus the spin corrections can all be 



obtained from the compact support pieces of the EW moments Eq. ( |2.19|) . We illustrate 
the procedure for computing the spin contributions by examining the 4-index EW multipole 
Eq. (pT9| c) 

/|$ = f r lj x k x l d 3 x. (Fl) 

J M 

Using Eq. fl2.9p and Eq. ( |2.5| ) we can write 

4S = f U V + (terms independent of velocity) + 0(^)1 X Wx . (F2) 

Terms which are independent of the fluid velocity will not contribute to the spin terms that 
we are computing here; they give non-spin terms which we have already calculated. Any 
spin terms that might result from the 0(pe 2 ) contributions will, in our convention, be at 
least 0{e 1 ^ 2 ) smaller, beyond the 2PN order at which we are working. We now write the 
source-point position and velocity as 

x i = x\ + x\ , (F3a) 
v* = v\ + v A , (F3b) 

where x % A is a suitably defined, PN-order, coordinate "center of mass" of body A and x\ 
is a coordinate displacement vector from the center of mass to the fluid element within the 
body. Similarly v A = dx A /dt is the coordinate velocity of the center of mass. (See e.g. 
f7H , |40|j41"|| for the definition of the center of mass.) 

Substituting Eq. ( |F3| ) into Eq. (|F2|) and integrating we obtain 

iSw = E m A vM + 2* iM& 4S7 , (F4) 

A 

where we have defined the spin vector by the formula 



L 



px\v\Sx = \e ki iS\ , (F5) 
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having assumed that J A px A v J 2d 3 x = (l/2)dP j (/dt = 0, where Pj is the body's intrinsic 
moment-of-inertia tensor. The first term in Eq. ( |F4| ) is the leading-order velocity-dependent 
term in Eq. (|4.26 ), and the second term is the spin-orbit correction to this multipole, of 
order e 1 / 2 smaller. In obtaining Eq. (F4) from Eq. ( |F2|) we have neglected a number of 
terms because (1) they vanish because of our assumption of spherical symmetry, (2) they 
have vanishing transverse-traceless projection, or (3) they are higher order in the bodies' 
small dimension (~ m), and therefore effectively of higher PN order. Such higher order 
moments can in principle be retained and incorporated into the framework. 

Keeping terms up to O(pe) in the source r*- 7 , and proceeding in precisely the same manner 
we can compute the spin-orbit contributions to the other EW multipoles 



TV 

1 EW 



^2 [^ax 1 a + x { X(y a x S A 
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(F6) 
(F7) 



Here again, the first terms are the leading-order non-spin contributions to the multipoles, 
Eqs. Q4.17D and ( [4.221) . The spin-orbit correction terms are respectively of order e 3//2 and 
e i/2 sma u er than the leading terms. 

In generating the expressions for the multipoles and waveforms, we must include spin 
corrections to the equations of motion. However, in the case of spinning bodies there is a 
delicate point to be considered in this procedure. The center of mass of body-A, denoted by 
xa used in our derivation of the multipole expressions turns out not to be precisely the same 
as the definition of the body's position used in the derivation of the conventional spin-orbit 
equations of motion, as given, say, by Damour [Q, or Eq. (|F14 ) below. The difference is 
related to the use of different so-called "spin supplementary conditions" which fix the center 
of mass of spinning bodies (see [4"D|,4~I] for a thorough discussion). We have previously shown 
f4~0|| that, to bring our center-of-mass definition into accord with that used in the equations 
of motion we need to shift the position of body-A in the following manner 
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2m A 



(v A x S A ) 



Performing this transformation replaces Eq. ( pq ) with 

,( l 7,r . S, C . V') 



1 EW 



[t^ax'a + 2x^(v A x S A ) 

A 



(F8) 



(F9) 



Since we are working only to 3/2PN order in the spin-orbit correction, the transformation 
Eq. ( |F8|) has no effect on the other multipoles. However if one were deriving the 2PN 
spin-orbit correction to the waveform [i.e. P 2 Q l g Q in Eq. (|6.10|) 1 it would be necessary to 
use the tranformation on Eq. ([FT]) as well. 
The spin pieces of Eqs. 



^9|), ([F7D and (jF4|) can just be added to their N-body point-mass 
counterparts in Section (|TV|), Eqs. ( |4. 1 T| j . ( |4. 22| ) and ( |4.26p , respectively. 

We now wish to restrict our attention to the two-body case and express our multipoles 
in terms of relative coordinates. The reduction parallels the 2-body (non-spin) reduction 
given in Section VI. We introduce the spin-corrections to the definition of the system center 
of mass, Eq. ( |4.16| c) (see p0| , |4l|| ), find the relation between the coordinates xx and x 2 and 
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the relative coordinate x corresponding to Eqs. (|6.2|) , and substitute into the two-body EW 
moments. It is useful to define two relative spin quantities 

Xs = U- 2 +^), (FlOa) 
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^=2U"4J- (F10b) 

With the spins normalized by the individual (masses) 2 , these vectors are essentially the 
vectorial sum and difference of the dimensionless angular-momentum (Kerr) parameters of 
the individual bodies. For orbital systems composed of two Kerr black holes or neutron 
stars these vectors will have a maximum magnitude of unity. Stability studies of rotating 
neutron stars show that the dimensionless angular momentum parameter is bounded above 
by 0.63 - 0.74, |72[ depending on the equation of state. Defining the vector spin quantities 



in this way also has the advantage that they are comparable in maximum magnitude to the 
other vectors that are used to form the terms in the waveform, namely n, N and v. As 
one computes the 2-body multipoles, the waveform, the energy flux, and the orbital phase 
evolution, the spins appear in many combinations with the masses. With the spin-quantity 
definitions as above, the reduced mass parameter 77 never appears in any denominators, so 



that the extreme mass ratio limit (77 — > 0) is always transparent in all expressions below [73 
This may seem like a minor aesthetic point, but it also means that the equations in the form 
we present them are suitable for stable numerical implementation with mass parameters free 
to roam from the equal mass case to the test mass case, and spin parameters free to roam 
independently of the mass choice from magnitude zero to unity. 



The spin corrections to the relation between x 1; x 2 and the relative coordinate x |4T 
take the form 

Xl = — x-mv x [x s (Sm/m) +%J , (Flla) 
m 

x 2 = x — mv x \x s (Sm/m) + x a ] ■ (Fllb) 

m 



Substituting these transformations into the leading order term in Eq. (p9|), we find that 
these spin-orbit corrections cancel, to the required order [compare Eq. (|6.3|)]. Substituting 
these definitions into the N-body multipoles gives the spin-orbit corrections to the two-body 
Epstein- Wagoner multipoles 

4w { so) = 4roV (v x Xs fx 3) , (F12a) 
4w { SO) = 2m VV )ifc [ (5m/m) Xs + *J , (F12b) 

I'iwiSO) = 4m 2 ??VV )m (V)xr • (F12c) 
These corrections can be added to the 2-body multipoles given in Section VI. STF multipoles 
can be projected from the EW multipoles using the formulae given in Appendix E. The 
results are 

r sTF { so) = f™V [2x*(v x Xs y - x X S ) J } STF , (F13a) 

3 r -i 
Jstf{so) = 2 m2 V [((<W m )Xs + XaY ^\ STF > (F13b) 



Jstf(so) = 4m V [sVtfJ STF ■ (F13c) 
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Eqs. ( |F13|) are in agreement with j|(],f|lj]. These spin-orbit contributions can be added 
to the STF multipoles given in Appendix E. It is interesting to note that the 4-index EW 
multipole I 1 ^ is needed to describe spin-dependence of the radiation, but there is no spin 
contribution from the 4-Index STF-multipole T^tf- The multipole 1%^ does contribute to 
the multipole TgTF an d Jstf through Eqs. (E2a) and (E2g). 

In order to derive the spin contributions to the waveform from the multipoles we must 
also augment the equations of motion [Eq. (|6.4p ] with spin-orbit and spin-spin contributions. 
These can be found in EDLET], and in our notation are given by 
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We now substitute our EW multipoles into Eq. (|2.18 ) and use the equations of motion to 
eliminate acceleration terms to obtain the final spin contributions to the waveform 
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Note that the spin-spin term comes entirely from the effects of the equations of motion. 
Thus we have computed the complete waveform, including leading-order spin effects, using 
our augmented EW formalism. The formalism can be extended to compute additional spin 
terms and other finite-size effects, such as the 2PN spin-orbit contribution to the waveform. 

Either by a direct computation starting with the waveform or by using the STF- 
multipoles in Eq. (|E5b|) we can compute the spin contributions to the rate of energy loss, 
Eq. (16121) , 
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Although they are not needed in our discussion, for completeness we include expressions 
for the precession of our spin vectors |40|J41[ 



mX s = n i x x s + n 2 x x a - 2(5m/m)x a x x s , 
m±a = n 2 xx„ + n 1 xx s - 2(1 - 2r]) Xs x x a . 



The precession vectors are given by 
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When spinning bodies are involved, the full gravitational- wave signal can become quite 
complicated; the orbital plane and the spin vectors of the individual bodies can precess, 
giving rise to a complicated modulation of the signal [41j,|74|. However in the special case 
when the spins are aligned (or anti-aligned) with the orbital angular momentum axis, the 
spin vectors and the orbital angular momentum vector do not precess [Eqs. ( F14j ), (|F18| ), 
( |F19| ) ] . In this special case there is a simple circular orbit solution to the equation of motion 
and it is straightforward to compute the spin contributions to the phase evolution. The spin 
contributions to orbital frequency can obtained from Eq. ( |F14| ), 
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where Xs and \ a now represent the projections of x s anc l X a onto the angular momentum 
axis. These quantities are positive when the spins are aligned in the same direction as the 
angular momentum axis and negative when they are anti-aligned. The orbital energy and 
energy flux take the simple form in the case of aligned spins and circular motion, 
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These spin corrections can be added to the non-spin formulae Eq. Q) and Eq. (Q. With 
these we can proceed as in Section VI to obtain the orbital angular velocity and orbital 
phase as explicit functions of time 
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Again, the spin-contributions can be inserted directly into Eqs. ( [7.11] ). (The definition of 
the dimensionless time T = rj(u/5m) is unchanged.) The explicit contributions to the + 
and x polarizations for this specialized circular orbit case can be obtained from Eq. ( [F15| ). 



In the notation of |48| they are given by 
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where x = mw and where the ". . ." represent the non-spin contributions given in [48 



In keeping with the notation used in [|S| the superscripts represent the post-Newtonian 
order and the physical nature of each term. The plus polarization spin-orbit and spin-spin 
contributions are 



H 



(1,50) 



H 



(3/2,50) 



H 



(2,55) 



-sini[(5m/m)x s + Xa] cos0 , 

- (1 + cos 2 i)[x s + {8m/m)xa\ + - 5cos 2 i)x s cos 20 

O 

-27/(1 +COS 2 0[( Xs ) 2 -(Xa) 21 



COS. 



and the cross polarization contributions are 
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We emphasize that these are only valid for quasi-circular orbits in the case where the 
the spins are aligned (or anti-aligned) with the orbital angular momentum vector. These 
restrictive assumptions about the configuration of the system suppress many of the intricate 
features of the waveform produced by spinning bodies [pH] , ^ . 

Figure 9 shows an inspiral waveform for the same system as in Figure 8 (1OM black hole 
and a 1.4M neutron star spiralling to coalescence), but in this case the objects are spinning. 
The spins are aligned with the orbital angular momentum axis. The spin contributions to 
both the waveform Eq. ( F23|) and the frequency evolution Eq. ( F22|) have been incorporated 
into the plot. The black hole has been given a spin of Sbh/^ 2 bh = 0-5 and the neutron 
star has S^s/ m2 Ns = 0-1 (*- e - Xs = 0.3 and Xa = 0.1). Notice the significant change in the 
frequency evolution; the system only sweeps to about 130 Hz in the same time it took for the 
non-spinning system to sweep to 180 Hz. Consequently, the peaks are not as closely bunched 
as they are in the non-spinning case. This slower orbital decay and frequency evolution is 
due to the dragging of inertial frames, which is inherent in the equations of motion and 
thus in our phase evolution equation Eq. ( |F22j ). At the left side of Figures 8 and 9, the 
waveforms are clearly in phase with each other, but after a few cycles they are out of phase. 
Since the phase evolution of the system is crucial in analysing gravitational waves from an 
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inspiral, it might seem that this sensitivity to spin in the phase evolution could be exploited 
and the spins of the bodies be determined with great accuracy. However, by leaving the 
spins the same but adjusting the masses slightly, we can recover the basic structure of the 
non-spinning case almost exactly. This is depicted in Fig. 10, in which the frequency sweep 
and the waveform itself are virtually identical to the non-spinning waveform in Fig. 8. This 



signal degeneracy in the spin and mass parameters has been previously noted in [24,25 



It is also interesting to notice that the inclusion of the spins virtually removes the jagged 
features from the troughs of the waves. 
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FIGURES 



FIG. 1. Past harmonic null cone C of the field point (i,x) intersects the near zone V in the 
hypersurface M. 

FIG. 2. Same as Fig. 1, for field point inside the near zone. 

FIG. 3. Taylor expansion of retarded time dependence on J\f results in multipole moments 
integrated over the spatial hypersurface M. 

FIG. 4. Two-dimensional hypersurfaces T formed by intersection of past null cone of field point 
with future null cones from the origin. Field point is in radiation zone. For u' from — oo to u — 21Z, 
F covers full 4ir solid angle around the origin. From u — 21Z to u, F terminates at boundary of the 
near zone M. 

FIG. 5. Same as Fig. 4, for field point in near zone. Integral over u' terminates at 
u' = u - 2K + 2r. 

FIG. 6. Fields contributing to A a @ at two representative points a and b on F have sources near 
same event u' at r = 0. Only orientation of near-zone source slice varies as angular integration 
moves around F. 

FIG. 7. Orientation of unit vectors defining + and x waveform polarizations. Direction of 
detector is N; p lies along line of nodes and is the origin for orbital phase angle 4>. 

FIG. 8. (a) Orbital frequency and (b) waveform for a 1.4M Q neutron star spiralling into a 1OM 
black hole plotted vs. time in seconds. Orbit is viewed edge-on, therefore only "+" -polarization is 
present. 

FIG. 9. Same configuration as Fig 8, but bodies are spinning. Both spins are aligned with 
orbital angular momentum axis. Angular momentum of black hole is = 0.5mf h and of neutron 
star is S ns = 0.5m^ s . Note frequency does not sweep as fast as non-spinning case because of 
dragging of inertial frames. 

FIG. 10. Spins are same as in Fig. 9, but heavier mass is now 12M Q . Frequency evolution is 
same as non-spinning case. Comparing this with Fig. 8 is an explicit demonstration of degeneracy 
in mass and spin parameters. 
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